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Abstract 

This  report  documents  the  CALGYP  model  which  is  designed  to  simulate  colcite 
and  gypsum  precipitation-dissolution  in  soils.  CALGYP  is  a  process  model  that 
is  easy  to  parameterize,  and  is  designed  for  long-term  simulations  (>  1000 
years).  The  CALGYP  model  has  five  components:  soil  parameterization, 
chemical  thermodynamic  relations,  deterministic  and  stochastic  rainfall  models, 
an  evapotranspiration  model,  and  subroutines  that  calculate  water,  calcium, 
and  sulfate  fluxes  through  the  soil.  The  stochastic  rainfall  model  is  based  on 
probability  distributions  for  interarrival  times  (days  between  rainfall  events)  and 
rainfall  amounts  and  is  designed  to  simulate  the  long-term  mean  annual  rainfall 
and  variability  in  annual  rainfall  for  specific  sites.  The  model  is  currently 
parameterized  for  seven  climatic  sites  in  the  desert  Southwest.  However,  climate 
(temperature  and  rainfall)  can  be  alfered  and  other  minerals  included,  which 
makes  the  CALGYP  model  potentially  applicable  across  a  wider  range  of 
environmental  conditions  including  freezing-thawing  systems.  A  separate 
program,  Rainmodule,  is  included  to  facilitate  inclusion  of  new  sites  and  to  alter 
rainfall  patterns  for  current  sites.  Instructions  for  utilization  and  a  FORTRAN-77 
source  code  listing  are  included  with  the  report. 


For  conversion  of  SI  metric  units  to  U.S./British  customary  units  of  measurement 
consult  ASTM  Standard  E380-89a,  Standard  Practice  for  Use  of  the  Irtternational 
System  of  Units,  published  by  the  American  Society  for  Testing  and  Materials, 
1916  Race  St.,  Philadelphia,  Pa.  19103. 
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CALGYP:  A  Simulation  Model  for  Calcite  and 
Gypsum  Precipitation-Dissolution  in  Soils 

GILES  M.  MARION 


INTRODUCTION 

Calcite  (CaCOa)  and  gypsum  (CaS04*2H20) 
are  common  secondary  minerals  precipitating  in 
both  hot-dry  and  cold-dry  soils  (Tedrow  1977, 
Harden et  al.  1991,  Marionet  al.  1993).  Soil  solution 
pH  in  calcareous  soils  is  largely  controlled  by  the 
solubility  of  calcite.  Quantifying  the  chemical  be¬ 
havior  of  semi-arid  to  arid  soils  in  both  hot  and 
cold  regions  requires  explicit  recognition  of  calcite 
and,  to  a  lesser  extent,  gypsum.  The  CALGYP 
model  was  originally  developed  to  predict  calcite 
and  gypsum  formation  in  desert  soils  of  the  South¬ 
west.  However,  climates  can  be  altered  and  other 
minerals  included  in  the  model,  which  makes  the 
CALGYP  model  potentially  applicable  across  a 
wider  range  of  environmental  conditions.  For  ex¬ 
ample,  CALGYP  could  be  used  to  assess  the  long¬ 
term  consequences  of  hazardous  waste  stabiliza¬ 
tion  in  a  CaCOs  matrix  (a  process  currently  under 
review  at  CRREL)  where  the  net  effect  over  time  is 
a  gradual  dissolution  and  removal  of  CaCOs.  Or 
CALGYP  could  be  used  to  simulate  salt  movement 
through  cold  regions  soils  such  as  those  that  exist 
along  the  Tanana  River  of  interior  Alaska,  which 
are  frequently  saturated  with  respect  to  both  cal¬ 
cite  and  gypsum  (Marion  et  al.  1993). 

Since  the  first  paper  based  on  this  model  (then 
called  CALDEP)  was  published  (Marion  et  al. 
1985),  copies  of  the  computer  code  have  been 
given  to  interested  scientists.  But  until  now,  this 
model  has  not  been  documented  in  a  publication 
that  would  render  the  model  more  readily  avail¬ 
able  to  interested  scientists.  The  objective  of  this 
report  is  to  document  CALGYP,  explaining  the 
theoretical  backgrormd,  the  practical  use,  the  limi¬ 


tations  of  the  model,  and  ho w  to  alter  the  model  for 
new  sites.  Also  included  are  the  FORTRAN  source 
code  listing;  a  copy  of  the  FORTRAN  program  on 
disk  is  available  from  the  author  on  request. 


MODEL  STRUCTURE 

Three  principles  guided  the  development  of  the 
CALGYP  model. 

1 .  The  model  must  be  process-based.  Only  at  the 
process  level  can  we  understand  fundamentally 
how  soil  processes  and  properties  interact  to  con¬ 
trol  soil  development. 

2.  Model  parameters  must  be  easily  estimated. 
This  would  facilitate  its  application  to  other  sites. 

3.  The  model  must  be  appropriate  for  long-term 
simulations  (>  1000  years).  This  latter  principle 
requires  simplification  to  improve  computational 
time. 

The  CALGYP  model  has  five  components:  soil 
parameterization,  chemical  thermodynamic  rela¬ 
tions,  stochastic  and  deterministic  rainfall  models, 
an  evapotranspiration  model,  and  subroutines 
which  calculate  water,  calcium,  and  sulfate  fluxes 
through  the  soil. 

Soil  parameterization 

CALGYP  is  a  compartment  model  that  can  be 
parameterized  to  contain  1  to  10  layers  (Fig.  1). 
Model  inputs  for  each  layer  include  layer  thick¬ 
ness,  bulk  density,  water  contents  initally  and  at 
soil  water  matric  suctions  of  0.01  MPa  (field  capac¬ 
ity)  and  1.5  MPa  (permanent  wilting  point),  initial 
soil  calcite  and  gypsum  contents,  initial  concentra¬ 
tions  of  soluble  Ca  and  SO4  -  S,  and  initial  soil  pH. 


Rainfall  Dust  Evapotranspiration 


Water,  Ca,  SO4 

Figure  1.  Schematic  diagram  of  the  CALGYP  model. 


A  model  for  soil  CO2,  developed  from  a  study  near 
Tucson,  Arizona  (Parada  et  al.  1983),  allows  CO2 
concentration  to  vary  spatially  and  seasonally 
(Marion  et  al.  1985).  Atmospheric  inputs  include 
the  Ca  and  SO4-S  content  of  rainfall  and  dust. 

Chemical  thermod5mamic  relations 

The  chemical  equilibrium  equations  used  in 
this  model  include  the  following: 


(C02)/Pc02  =  Ki 

(1) 

pKi  =  1.14  + 0.0131  f 

(2) 

(H+)(HC03-)/(H20)(C02)  =  K2 

(3) 

pK2  =  6.54  -  0.0071 1 

(4) 

(H+)(C032-)/(HC03-)  =  K3 

(5) 

pK3  =  10.59-0.0102  t 

(6) 

(Ca2+)(C032-)  =  K4 

(7) 

pK4  =  7.95  +  0.0125  t 

(8) 

(Ca2+)(S042-)(H20)2  =  K5 

(9) 

pKs  =  4.62  +  0.0006 1 

(10) 

(Ca2+)(S042-)/(CaS04°)  =  Ke 

(11) 

pKe  =  2.23  +  0.0019 1 

(12) 

where  pK  is  the  negative  logarithm  of  the  equilib- 
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riirai  constant,  t  is  temperature  (°C),  and  parenthe¬ 
ses  refer  to  activities.  The  equilibrium  constants 
and  their  temperature  dependencies  were  esti¬ 
mated  from  equilibrium  data  over  the  tempera¬ 
ture  range,  0  to  40  °C  (Garrels  and  Christ  1965) .  The 
intercept  term  in  eq  8  was  selected  to  yield  a  pK  of 
8.26  at  25  °C,  which  was  the  mean  of  several  Bk2 
horizon  samples  from  the  desert  LTER  site  equili¬ 
brated  at  25  °C  tmder  a  fixed  CO2  concentration 
(500  ppm)  for  10  days  (Marion  et  al.  1990). 

Single-ion  activities  (a)  of  ions  having  aqueous 
solution  concentrations  of  c  are  estimated  by 

a  =  yc  (13) 

where  yis  the  single-ion  activity  coefficient,  which 
was  estimated  with  the  Davies  equation  (Davies 
1962) 

log  Y  =  -  A  z2  [Vf/(1.0  +  V/)  -  0.3 1]  (14) 

where  z  is  the  ionic  valence,  I  is  the  ionic  strength 
which  is  defined  as 

J  =  0.5E(CiZi2)  (15) 

and  A  is  the  Debye-Hiickel  constant  which  is  given 
by 


A  =  0.4918  +  6.6098  x  10-4  t  +  5.0231  x  10-6  ^2 

(16) 

over  the  temperature  range:  0-40  °C  (Robinson 
and  Stokes  1965). 

Critical  to  the  performance  of  these  models  is 
the  mechanism  for  estimating  soil  solution  pH, 
because  of  the  strong  influence  of  pH  on  CaCOa 
solubility.  For  a  pure  CaC03-CaS04  solution  in 
the  pH  range  of  7-8.5,  the  following  charge  bal¬ 
ance  holds: 

2[Ca2+]  =  [HCO3-]  +  2[C032-]  +  2[S042-] 

(17) 

where  brackets  refer  to  concentrations.  For  a  sys¬ 
tem  in  equilibrium  with  solid  CaC03,  this  equa¬ 
tion  can  be  rewritten  as: 

2K4(H+)^  _  K1K2PC02  , 

YCaK3K2KlI’C02  (H‘^)YhC03 


2Ki  K2K3PC02 
(H+)^  YC03 


+  2  [SO42-]  . 


Given  the  Pco2  (partial  pressure  of  CO2)  and  the 
SO4-S  concentration,  eq  1 8  is  solved  for  (H-^),  which 
is  then  used  to  control  CaC03  solubility. 

For  a  system  in  equilibrium  with  both  solid 
calcite  and  gypsum,  the  sulfate  term  in  brackets 
(eq  17  and  18)  may  be  replaced  by 


K5K1K2K3PC02  .9. 

K4YS04(Hf 

Compared  to  equilibrium  with  pure  calcite,  simul¬ 
taneous  equilibrium  of  gypsum  and  calcite  signifi¬ 
cantly  depresses  the  equilibrium  pH  and  increases 
the  equilibrium  Ca  concentrations  (Fig.  2).  In  the 
CALGYP  model,  eq  18  is  used  to  estimate  solution 
pH  (after  CaC03  begins  to  precipitate),  so  that  soil 


(18) 
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Figure  2.  Theoretical  (A)  pH  and  (B)  Ca  concentrations 
for  pure  calcite  and  calcite-gypsum  solutions  as  func¬ 
tions  of  the  partial  pressure  of  carbon  dioxide  at  25  °C. 


Precipitation  (cm/day) 


Figure  3.  Cumulative  probabilities  for  interarrival  daps  and  daily  rainfall  for  the  Tucson  site. 


solution  pH  is  a  continuous  function  of  soil  solu¬ 
tion  sulfate  concentration.  Empirical  justification 
for  this  assumption  was  demonstrated  in  a  previ¬ 
ous  paper  (Marion  and  Schlesinger  in  press). 

If  soil  solution  sulfate  concentrations  are  "zero/' 
CALGYP  will  by-pass  the  sulfate  equilibrium  rou¬ 
tines  and  will  calculate  chemical  equilibrium  for  a 
pure  calcite  system.  One  can  also  effectively  by¬ 
pass  calcite  equilibrium  and  its  pH  dependence 
(eq  18)  by  assigning  the  initial  soil  an  arbitrarily 
low  pH  of  5 .0.  At  this  low  pH,  it  is  doubtful  that  Ca 
concentrations  will  ever  build  up  to  the  point 
where  calcite  would  precipitate  given  normal  soil 
CO2  concentrations. 

Rainfall  models 

The  stochastic  rainfall  model  controls  input  of 
water  and  is  based  on  probability  distributions  for 
interarrival  times  (the  number  of  days  between 
rainfall  events)  and  the  rainfall  amounts  for  spe¬ 
cific  seasons  at  specific  sites  (Fig.  3) .  Sites  currently 
in  the  model  include  Yuma,  Phoenix,  and  Tucson 
in  Arizona;  Albuquerque,  Roswell,  and  Clayton  in 
New  Mexico;  and  El  Paso,  Texas.  A  random-num¬ 
ber  generator  is  used  to  select  the  interarrival 
times  and  rainfall  amounts  for  each  year  from  the 
cumulative  probability  distributions.  This  stochas¬ 
tic  model  is  designed  to  reproduce  the  long-term 
average  annual  rainfall  and  the  variability  in  an¬ 
nual  rainfall  for  a  specific  site  (Marion  et  al.  1985). 

In  addition  to  the  stochastic  rainfall  model, 
CALGYP  also  includes  an  option  for  a  determinis¬ 
tic  rainfall  model  that  allows  the  user  to  specify  the 
yearly  rain  dates  and  rainfall  amounts  for  a  spe¬ 


cific  site.  The  same  rainfall  pattern  is  then  used 
year-after-year. 

Evapotranspiration  model 

The  evapotranspiration  model,  which  is  prima¬ 
rily  a  function  of  temperature,  controls  the  loss  of 
water  and  consists  of  three  steps.  First,  potential 
evapotranspiration  is  calculated  using  Thom- 
thwaite's  equation  (Thomthwaite  1948,  Marion  et 
al.  1985).  Second,  Thomthwaite's  potential  evapo¬ 
transpiration  is  converted  to  pan  evaporation  us¬ 
ing  a  derived,  empirical  relationship  with  tem¬ 
perature  for  Southwestern  deserts  (Fig.  4).  And 
third,  actual  evapotranspiration  is  calculated  as  a 
function  of  soil  moisture  and  pan  evaporation 
between  field  capacity  (0.01  MPa)  and  permanent 
wilting  point  (1.5  MPa)  (Fig.  5).  Calibration  of  the 
third  step  is  based  on  field  measurements  from  a 
Larrea  tridentata  (creosote  bush)  site  at  the  Jornada 
Desert  Long-Term  Ecological  Research  site  near 
Las  Cruces,  New  Mexico  (Marion  et  al.  1985). 
Water  loss  is  assumed  to  occur  at  the  potential  rate 
(ratio  =  1.00)  in  the  upper  45.4  %  of  the  available 
moisture  range;  in  the  lower  54.6  %  of  the  range, 
water  loss  is  a  linear  function  of  soil  moisture  (Fig. 
5). 

CALGYP  includes  options  to  change  the  cli¬ 
matic  variables  of  rainfall  and  temperature.  For  a 
temperature  change  for  the  Southwestern  desert 
sites,  CALGYP  assumes  that  monthly  pan  evapo¬ 
ration  is  decreased  (increased)  by  an  amount  pro¬ 
portional  to  the  mean  annual  temperature  de¬ 
crease  (increase)  (Fig.  6).  For  a  fuller  discussion  on 
altering  climate  within  and  outside  the  desert 
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Figure  4.  Pan  evaporation/Thomthwaite  potential  evapotranspiration  ratio  as  a  func¬ 
tion  of  mean  monthly  temperature  for  Southwestern  desert  sites. 


Figure  5.  Relationship  of  the  actual  evapotranspiration/pan  evaporation  ratio  as  a  function  of  soil 
moisture  content. 
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Figure  6.  Annual  pan  evaporation  (calculated)  as  a  function  of  mean  annual 
temperature. 


Southwest,  see  Altered  Climate  and  Chemistry. 

Water  and  solute  flux 

A  daily  time-step  is  used  to  assess  the  fluxes  of 
water  and  solutes  through  the  soil.  All  rainfall  is 
assumed  to  enter  the  uppermost  soil  layer  (Fig.  1); 
the  model  ignores  vegetative  interception  of  rain¬ 
fall  and  surface  runoff.  Only  saturated  flow  of 
water  is  considered  in  these  models.  If  the  water¬ 
holding  capacity  of  a  layer  is  exceeded  (i.e.,  water 
content  >  field  capacity),  excess  water  moves  into 
progressively  deeper  layers.  Water  flux  beyond 
the  base  of  the  soil  profile  is  treated  as  leachate  and 
is  assiuned  lost  from  the  system.  Solutes  are  as¬ 
sumed  to  move  with  the  mass  flow  of  watier.  Water 
that  enters  a  given  layer  is  mixed  with  the  pre¬ 
existing  water  and  solutes  are  equilibrated  chemi¬ 
cally  with  the  solid  and  gas  phases.  Therefore,  the 


excess  water  that  passes  through  a  given  layer 
contains  an  equilibrated  concentration  of  solutes 
before  passing  to  the  deeper  layers.  During  drying 
cycles,  water  is  first  extracted  from  the  surface 
layer  and  then  from  progressively  deeper  layers 
using  the  evapotranspiration  model  previously 
mentioned. 


C  ALG  YP  FLOWCHART 

After  initialization  of  soil  properties  and  calcu¬ 
lation  of  monthly  potential  evaporation,  CALGYP 
cycles  along  three  time  steps  (Fig.  7).  Water  loss  is 
calculated  on  a  daily  time  step.  Whenever  a  rain 
event  occurs,  CALGYP  cycles  through  chemical 
equilibrium  and  water  flux  routines  to  calculate 
the  fluxes  of  water  and  solutes  through  the  soil 
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Figure  7.  CALGYP  program  flowchart. 


profile.  The  chemical  equilibrium  routines  are 
addressed  twice  for  each  rain  event.  First,  calcite 
and  gypsum  precipitation  are  calculated  over  the 
previous  drying  cycle  up  to  the  rain  event  day. 
Then,  dissolution  and  flow  of  calcite  and  gypsum 
are  calculated  after  the  water  flux  determination. 
An  armual  cycle  includes  calculation  of  the  sto¬ 
chastic  rain  pattern,  the  daily  and  rain  event  cycles, 
and  print  statements. 

THE  FORTRAN  PROGRAM 

The  original  version  of  the  CALGYP  program 
was  called  CALDEP  because  it  dealt  exclusively 
with  calcite  deposition  and  was  written  in  HP  Basic 
(Marion  et  al.  1 985) .  This  version  required  about  45 


seconds  of  computer  time  to  simulate  a  year  of  soil 
development  on  an  HP  9816  microcomputer.  This 
program  was  translated  into  a  TrueBASIC  version 
for  operation  on  an  Apple  Macintosh  Ilex  com¬ 
puter.  When  gypsum  was  added  to  the  model,  the 
name  was  changed  to  CALGYP.  The  CALGYP 
program  required  11  to  22  seconds/year  on  16- 
MHz  and  25-MHz  Apple  Macintosh  computers, 
respectively.  Recently  CALGYP  was  translated  to 
FORTRAN,  primarily  to  increase  portability  to 
other  computer  systems  and  to  improve  computa¬ 
tional  time.  The  present  FORTRAN  version  re¬ 
quires  0.22  seconds/year  running  on  a  33-MHz 
Apple  Macintosh  Quadra  800 using  the  MacFortran 
II  compiler  optimized  for  the  Motorola  68040  pro¬ 
cessor  (Absoft  1993).  To  simulate  10,000  years  of 
change  now  requires  37  minutes  with  the  FOR- 
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TRAN  version  on  a  fast  microcomputer  compared 
to  5.2  days  with  the  original  HP  Basic  version  on  a 
slower  microcomputer.  This  200-fold  improve¬ 
ment  in  computational  speed  is  a  major  enhance¬ 
ment  because  this  program  must  be  able  to  simu¬ 
late  slow  change  over  very  long  periods. 

A  major  effort  was  made  to  adhere  to  the  FOR¬ 
TRAN  77  ANSI  Standard  in  the  translation.  For 
example,  CALGYP  uses  an  80-column  statement 
line  where  columns  1-5  are  statement  labels,  col¬ 
umn  6  is  a  continuation  column,  coluirms  7-72  are 
used  for  statements,  and  columns  73-80  are  used 
for  commentary.  To  facilitate  the  translation,  vari¬ 
able  names  were  retained  between  programs 
except  to  shorten  names  with  >  six  characters  in 
the  TrueBASIC  version  to  <  six  characters  in  the 
FORTRAN  version.  FORTRAN  assumes,  rmless 
otherwise  declared,  that  variables  beginning  with 
the  letters  "I,  J,  K,  L,  M,  and  N"  are  integer  vari¬ 
ables  and  all  other  initial  letters  refer  to  real  vari¬ 
ables.  In  order  to  retain  the  same  variable  names 
between  programs,  this  restriction  required  for¬ 
mal  declaration  of  some  real  and  integer  variables 
in  the  FORTRAN  version.  Also,  in  a  few  cases, 
integer  arithmetic  had  to  be  rewritten  in  terms  of 
real  variables  to  preserve  meaning  between  pro¬ 
grams. 

The  CALGYP  program  consists  of  a  main  pro¬ 
gram,  also  called  CALGYP,  and  eight  subroutines, 
four  for  climate  and  four  for  chemical  equilibrium. 
In  addition,  a  separate  program  called  Rainmodule 
is  included  in  order  to  facilitate  altering  rainfall 
patterns.  CALGYP  and  Rainmodule  use  one  exter¬ 
nal  library  routine  called  "rand"  (UNIX  library) 
which  selects  random  numbers  between  0.0  and 
1.0.  Program  listings  are  included  in  Appendices 
A  and  B. 

In  contrast  to  the  "modular"  approach  where 
the  main  program  primarily  "calls"  other  subpro¬ 
grams  (Nyhoff  and  Leestma  1985,  Plummer  et  al. 
1988),  the  CALGYP  main  program  directly  per¬ 
forms  many  important  tasks.  CALGYP:  1)  controls 
screen  queries,  2)  contains  the  bulk  of  the  data 
used  by  the  program,  3)  initiates  variables,  4) 
calculates  moisture  flow,  5)  mixes  solutions,  6) 
updates  soil  properties,  and  7)  prints  results. 

The  four  climate  subroutines  are  Detrain, 
Raindate,  Rainfall,  and  Seasons.  Detrain  allows  the 
user  to  specify  rain  dates  and  rainfall  amounts  for 
a  given  year  at  a  specific  site.  This  yearly  pattern  is 
then  reused  year  after  year  (the  "deterministic" 
rain  model).  The  subroutines  Raindate,  Rainfall, 
and  Seasons  are  used  to  develop  the  "stochastic" 
rain  model,  which  varies  rain  dates  and  rainfall 


amormts  over  annual  cycles  based  on  probability 
distributions  and  a  random-number  generator. 
Raindate  and  Rainfall  calculate  the  armual  rain 
dates  and  rainfall  amoimts  for  a  specific  site,  re¬ 
spectively.  Seasons  calculates  the  proper  season 
(e.g.,  winter,  spring,  summer)  for  a  given  day  of 
the  year  and  site.  This  controls  the  proper  prob¬ 
ability  distributions  used  in  the  Raindate  and  Rain¬ 
fall  subroutines. 

The  four  chemical  subroutines  are  Constants, 
Ion,  Hact,  and  Equil.  Constants  calculates  the  equi¬ 
librium  constants  (eq  1-12)  and  the  Debye-Hiickel 
constant  "A"  (eq  14  and  16)  for  a  specific  tempera- 
txue.  Mean  monthly  air  temperatures  for  a  specific 
site  are  used  for  these  temperature  calculations. 
Ion  calculates  the  single-ion  activity  coefficients  of 
xmivalent  and  bivalent  ions  using  the  Davies  equa¬ 
tion  (eq  14).  Hact  calculates  the  H  activity  using  eq 
18.  Equil  calculates  the  equilibrium  composition  of 
solutions  accounting  for  the  CaS04°  ion-pair  and 
the  precipitation  and  dissolution  of  calcite  and 
gypsum. 


PROGRAM  INPUT 
CALGYP  program 

CALGYP  input  occurs  primarily  through  screen 
queries,  DATA  statements,  and  two  subroutines 
{Detrain  and  Constants).  CALGYP  is  designed  to 
work  interactively  with  on-screen  prompts  re¬ 
questing  information  and  options  to  control  the 
simulation  (SeeTable  1).  Theinformationrequested 
includes: 

(1)  Site  Selection?  This  prompt  requests  that 
the  user  select  one  of  seven  sites  by  entering 
the  site  designation  (Table  1),  then  the  RE¬ 
TURN  key. 

(2)  Current  Climate  or  Altered  Climate?  En¬ 
ter  proper  designation,  then  RETURN.  If 
Altered  Climate  is  selected,  then: 

(a)  Change  in  temperature  (°C).  Enter  tem¬ 
perature  change  (+  for  increase,  -  for  de¬ 
crease). 

(b)  Fractional  change  in  rainfall  amount 
during  drier  climate.  If  rainfall  is  to  de¬ 
crease  10%,  enter  0.10,  then  RETURN. 
Wetter  climates  require  altering  the  prob¬ 
ability  distributions  (See  Altered  Climate 
and  Chemistry). 

(3)  Title?  Enter  any  alphanumeric  title  up  to 
50  characters,  then  RETURN. 

(4)  Number  of  Soil  Horizons?  CALGYP  can 
work  with  1  to  10  horizons  (layers).  This 


8 


Table  1.  An  example  of  the  CALGYP  screen  query  input. 


Select  Site: 

1  =  El  Paso,  TX  (21.6cm) 

2  =  Albuq.,  NM  (21.1cm) 

3  =  Clayton,  NM  (37.8cm) 

4  =  Roswell,  NM  (31.6cm) 

5  =  Yuma,  AZ  (8.5cm) 

6  =  Phoenix,  AZ  (18.9cm) 

7  =  Tucson,  AZ  (28.4cm) 

7 

Climate  Option?  Current  Climate  =1,  Altered  Climate  =  2 
2 

Enter  Delta  Temperature=(®C,default=0.0,no  change) 

-5 

Enter  Fractional  Change  in  Rainfall  Amount  during 
Drier  Climate  =  (default  =0.0,  no  change) 

0.1 

Title? 

TUCSON  SIMULATION  OF  ALTERED  CLIMATE 
Number  of  Soil  Horizons  (Max=10) ? 

10 

Years  to  Simulate? 

1 

Print  Interval? 

1 

Deterministic  (1)  or  Stochastic  (2)  Rainfall  Model? 

2 

Seed  for  Random  Number  Generator? 

123 


number  must  not  exceed  the  dimension  of 
the  soU  horizon  properties  specified  in  the 
DATA  statements. 

(5)  Years  to  Simulate?  Enter  the  number  of 
years  that  you  wish  this  simulation  to  nm, 
then  RETURN. 

(6)  Print  Interval?  At  what  yearly  interval  do 
you  wish  intermediate  results  printed.  For 
example,  for  a  1000-year  simulation,  you 
may  wish  to  print  proffle  descriptions  at  200- 
year  intervals. 

(7)  Deterministic  or  Stochastic  Rainfall 
Model?  Enter  designation  for  appropriate 
rainfall  model.  If  stochastic  is  selected,  then: 

fal  Seed  for  Random  Number  Generator? 
Enter  integer  of  1-5  digits.  Using  the  same 
Seed  in  nms  will  produce  the  same  se¬ 
quence  of  random  numbers,  which  is  par¬ 
ticularly  useful  during  the  development 
and  debugging  stages. 

Table  1  is  an  example  of  a  screen  query  that 
selects  the  Tucson  site  with  altered  climate  (5  °C 
colder  with  10%  less  rainfall).  The  number  of  soil 
horizons  is  10,  and  the  simulation  nms  for  one  year 
with  a  printout  at  one  year.  The  stochastic  rainfall 
model  was  selected  with  a  Seed  of  123  for  the 
random  number  generator. 

Most  of  the  data  used  to  parameterize  the  model 
are  stored  in  DATA  statements  within  the  CALGYP 


program.  These  DATA  statements  start  on  p.  2  of 
the  program  listing  (App.  A).  To  make  changes  in 
these  data  requires  entering  and  altering  the  source 
code. 

The  first  DATA  statement  lists  the  mean  monthly 
air  temperatures  (°C)  for  the  seven  sites.  Commen¬ 
tary  in  columns  73-80  identifies  the  specific  site. 
The  temperature  data  are  followed  by  several 
DATA  statements  that  specify  the  parameters  for 
the  stochastic  rainfall  model  by  site  and  season. 
Five  of  the  sites  have  two  rainfall  seasons  (winter 
and  summer).  Phoenix  and  Tucson  have  three 
rainfall  seasons  (winter,  spring,  summer).  Details 
on  the  derivation  of  these  distributions  were  dis¬ 
cussed  in  Marion  et  al.  (1985) .  The  variables  D(I,J,K) 
and  R(I,J,K)  are  probabilities  and  corresponding 
daily  rainfall  amoimts.  For  example,  the  probabil¬ 
ity  is  0.837  of  a  daily  rainfall  amount  <  1.00  cm  for 
El  Paso  in  the  winter  (App.  A).  For  this  site  and 
season,  a  random  number  of  0.600  would  fall 
between  rainfall  amounts  of  0.25  and  0.50  cm; 
linearly  interpolating  gives  0.29  cm  of  rain  at  a 
probability  of  0.600.  In  a  few  cases,  dummy  prob¬ 
abilities  (>  1.00)  are  present  in  order  to  complete 
matrices.  The  only  physically  meaningful  prob¬ 
abilities  are  <  1 .00.  For  example,  the  only  meaning¬ 
ful  probabilities  [D(I,J,K)]  for  Phoenix  in  the  spring 
are  0.000  to  1 .000 ;  the  values  ranging  from  1 . 1  to  1 .6 
are  dummy  variables  as  are  the  corresponding 
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rainfall  amounts  ranging  from  1 .25  to  3.07  cm.  The 
variables  Freq  (I  J,K)  and  Inter  (I,J,K)  are  the  prob¬ 
abilities  and  corresponding  interarrival  days  (the 
number  of  days  between  rainfall  events).  For  El 
Paso  in  the  winter,  the  probability  is  0.765  of  <  10 
days  between  rain  events.  For  this  site  and  season, 
a  random  number  of  0.600  corresponds  to  an 
interarrival  time  of  5.8  days. 

Subsequent  DATA  statements  define  the  basic 
soil  properties  of  the  system.  Properties  that  must 
be  specified  by  horizon  include  thickness,  bulk 
density,  soil  water  concentrations  at  0.01  MPa 
(field  capacity),  1 .5  MPa  (permanent  wilting  point), 
and  at  first,  initial  CaCOa  and  soluble  Ca,  initial 
soil  solution  H  activity  and  soil  atmospheric  CO2 
concentration,  and  initial  CaS04  and  soluble  SO4. 
The  program  allows  the  user  to  change  monthly 
CO2  concentrations  in  the  soil  profile  by  using  a 
single  monthly  multiplier  for  the  entire  profile.  In 
the  present  program  ( App .  A),  the  winter  months' 
(November-February)  CO2  concentrations  are  in¬ 
creased  by  58.8%  during  the  other  months  (Marion 
et  al.  1985). 

The  final  data  statements  define  the  atmospheric 
properties  that  include  the  dust  Ca  and  SO4  con¬ 
tents  and  the  rain  Ca  and  SO4  concentrations. 
CALGYP  only  considers  precipitation  and  disso¬ 
lution  of  calcite  and  gypsum.  Other  minerals  that 
might  contribute  to  Ca  and  SO4  input  are  ignored, 
which  is  a  reasonable  assumption  for  Southwest¬ 
ern  desert  soils  because  of  the  generally  low  rates 
of  mineral  weathering  (Gile  et  al.  1981,  Machette 
1985,  Harden  et  al.  1991,  McFadden  et  al.  1991).  If 
weathering  of  minerals  other  than  calcite  and  gyp¬ 
sum  were  important  and  these  rates  were  known, 
they  could  be  included  with  dust  inputs. 

The  subroutine  Detrain  contains  ^e  rain  dates 
and  rainfall  amounts  for  the  deterministic  rain 
model.  This  subroutine  requires  specification  of 
the  rain  dates  (Julian  day)  and  rainfall  amounts 
(cm)  in  the  appropriate  DATA  statements  {Raindy 
and  Raina).  Also,  the  Dimension  statement  and  the 
variable,  Sumfre,  must  be  changed  to  reflect  the 
total  number  of  rain  events/ year.  Constants  con¬ 
tains  the  chemical  thermodynamic  constants  and 
their  temperature  dependence. 

Rainmodule  Program 

The  Rainmodule  program  was  designed  to  fa¬ 
cilitate  inclusion  of  new  sites  with  new  rainfall 
patterns  and  for  altering  rainfall  for  current  sites. 
The  structure  of  this  program  is  similar  to  that 
used  for  the  stochastic  model  in  CALGYP.  Rim- 


ning  this  model  "as  is"  will  result  in  rainfall  pat¬ 
terns  identical  to  those  used  in  CALGYP.  By  trial 
and  error,  one  can  alter  either  the  frequency  distri¬ 
butions  or  the  variables  "Newj"  or  "Numb"  to 
increase  or  decrease  rain  dates  or  rainfall  amoimts. 
In  general,  the  frequency  distributions  are  derived 
from  short-term  records  (e.g.,  3-10  years)  which 
may  be  wetter  or  drier  than  normal.  The  variables 
Newj  and  Numb  are  used  to  increase  or  decrease 
the  number  of  rain  dates  per  year,  respectively,  in 
order  to  match  the  long-term  (40-100  year  climate 
record)  mean  armual  rainfall  for  a  specific  site. 
Newj  adds  rain  dates  increasing  rainfall;  Numb 
assigns  0.0  to  rainfall  amounts,  effectively  elimi¬ 
nating  dates  and  decreasing  rainfall.  Manipula¬ 
tion  of  the  frequency  distributions  will  be  dis¬ 
cussed  under  Altered  Climate  and  Chemistry. 


PROGRAM  OUTPUT 
CALGYP  program 

Output  from  CALGYP  is  printed  to  a  file  called 
"CaData"  in  three  major  blocks.  Table  2  contains 
an  example  for  a  1000-year  simulation  with  inter¬ 
mediate  results  printed  at  200-year  intervals. 

The  first  block  contains  the  title  and  a  few  key 
run  options  followed  by  a  climate  summary 
(monthly  temperatures  and  calculated  potential 
evapotranspiration),  the  initial  soil  profile,  and 
atmospheric  chemical  drivers. 

The  second  block,  repeated  several  times,  con¬ 
tains  intermediate  soil  profile  properties  at  the 
"Print"  interval.  These  properties  include  CaCOs 
and  CaS04  contents,  the  bulk  density  (which  in¬ 
creases  as  mineral  salts  precipitate),  the  total 
amoimt  of  water,  Ca,  and  S  leached  past  the  base 
of  the  soil  profile  through  the  specified  year,  and 
the  annual  rainfall  and  evapotranspiration  for  the 
specified  year. 

The  third  block  includes  a  few  additional  soil 
profile  properties  at  the  end  of  the  simulation, 
such  as  moisture  content,  H  ion  activity,  and  soluble 
Ca  and  SO4.  Also  total  rain  and  total  evapotranspi¬ 
ration  for  the  entire  simulation  are  printed.  And 
finally  residual  dust  Ca  and  SO4  are  printed.  The 
model  accumulates  dust  at  the  soil  surface  on  a 
daily  basis  and  washes  it  into  the  profile  whenever 
it  rains.  "Residual  dust"  refers  to  the  amormt  of  Ca 
and  SO4  that  has  accumulated  as  dust  at  the  soil 
surface  between  the  last  rain  and  the  end  of  the 
year  for  the  last  year  of  the  simulation.  This  final 
data  block  may  be  useful  intrinsically  (e.g.,  for 
mass  balance  calculations)  as  well  as  for  providing 
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Table  2.  Sample  output  of  the  CALGYP  model  for  a  1000-year  simulation  for  the  Tucson  site. 


4/18/94  10:37 


Macintosh  HD:MPW:CaData 


TUCSON  1000  YR  SIMULATION 
SITE  =  7 

Stochastic  Rainfall  Model  with  Seed  =  123 


CLIMATIC  SUMMARY 


Month 

Temp (C) 

1 

10.2 

2 

11.8 

3 

14.3 

4 

18.1 

5 

22.6 

6 

27.8 

7 

30.0 

8 

28.9 

9 

26.7 

10 

20.8 

11 

14.6 

12 

10.8 

Pot.Evap. (g/cm2) 
10.90 
11.77 
16.26 
20.44 
27.47 
35.39 

41.19 

28.19 
24.37 
■18.54 
12.34 

9.09 


INITIAL  SOIL  PROFILE 


Hori 

Thick 

BD 

C02 

H 

CaC03 

Ca 

CaS04 

S04 

cm 

g/cm3 

atm 

gCa/cm2 

gCa/cm2 

gS/cm2 

gS/cm2 

1 

10.0 

1.44 

0.379E-03 

O.lOOE-07 

O.OOOE-fOO 

0.864E-05 

O.OOOE+00 

0.144E-04 

2 

10.0 

1.44 

0.687E-03 

O.lOOE-07 

O.OOOE+00 

0.864E-05 

O.OOOE+00 

0.144E-04 

3 

10.0 

1.44 

0.976E-03 

O.lOOE-07 

O.OOOE+00 

0.864E-05 

O.OOOE+00 

0.144E-04 

4 

10.0 

1.44 

0.128E-02 

O.lOOE-07 

O.OOOE-t-00 

0.864E-05 

O.OOOE+00 

0.144E-04 

5 

10.0 

1.44 

0.160E-02 

O.lOOE-07 

O.OOOE-fOO 

0.864E-05 

O.OOOE+00 

0.144E-04 

6 

10.0 

1.44 

0.191E-02 

O.lOOE-07 

O.OOOE+00 

0.864E-05 

O.OOOE+00 

0.144E-04 

7 

10.0 

1.44 

0.214E-02 

O.lOOE-07 

O.OOOE-i-00 

0.864E-05 

O.OOOE+00 

0.144E-04 

8 

10.0 

1.44 

0.227E-02 

O.lOOE-07 

O.OOOE+00 

0.864E-05 

O.OOOE+00 

0.144E-04 

9 

10.0 

1.44 

0.241E-02 

O.lOOE-07 

O.OOOE+00 

0.864E-05 

O.OOOE+00 

0.144E-04 

10 

10.0 

1.44 

0.254E-02 

O.lOOE-07 

O.OOOE+00 

0.864E-05 

O.OOOE+00 

0.144E-04 

ATMOSPHERIC  CHEMICAL  CONDITIONS 

Dustca <gCa/cm2/yr)  Precca (mgCa/1)  Dusts (gS/cin2/yr)  Precs(mgS/l) 
0.8614E-04  .00  0.7993E-05  .00 


YEAR  =  200 

CAC03(gCa/cni2)  =  O.OOOOOE+00  0.16325E-04  0.13443E-03  0.55263E-02  0.42430E-02 
CAC03(gCa/cin2)  =  0.24938E-02  0.13452E-02  0.45815E-03  0.23268E-03  0.17588E-03 

CAS04  (gS/cm2)  =  0.00000E■^00  O.OOOOOE-fOO  O.OOOOOE+00  O.OOOOOE+00  O.OOOOOE+00 
CAS04 (gS/cm2)  =  O.OOOOOE+00  O.OOOOOE+00  0.60269E-04  0.19028E-03  0.30192E-03 

Bulk  Density (g/cm3)  =  1.4400  1.4400  1.4400  1.4414  1.4411 

Bulk  Density (g/cni3)  =  1.4406  1.4403  1.4401  1.4402  1.4402 

Leach(cm)  =  0.6681E+01  Leacca (gCa/cm2)  =0.47531E-03  Leachs (gS/cm2)  =0.17596E-03 
Annual  Rain (cm)  =  27.899  Annual  Evap(cm)  =  27.858 


YEAR  =  400 

CAC03(gCa/cm2)  =  O.OOOOOE+00  0.63449E-04  0.65553E-03  0.10642E-01  0.89438E-02 
CAC03(gCa/cm2)  =  0.52801E-02  0.28194E-02  0.94581E-03  0.34664E-03  0.17595E-03 

CAS04  (gS/cm2)  =  O.OOOOOE+00  0 .  OOOOOE-t-OO  0 .  OOOOOE-t-00  O.OOOOOE+00  O.OOOOOE+00 
CAS04 (gS/cm2)  =  O.OOOOOE+00  O.OOOOOE+00  0.38755E-03  0.12687E-02  0.30205E-03 
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Table  2  (Cont'd).  Sample  output  of  the  CALGYP  model  for  a  1000-year  simulation  for  the 
Tucson  site. 


4/18/94  10:37 


Macintosh  HD:MPW:CaData 


Bulk  Density (g/cm3)  =  1.4400  1.4400  1.4402  1.4427  1.4422 

Bulk  Density (g/cm3)  =  1.4413  1.4407  1.4404  1.4408  1.4402 

Leach (cm)  =  0.6681E+01  Leacca (gCa/cm2)  =0.47531E-03  Leachs (gS/cm2)  =0.n596E-03 
Annual  Rain (cm)  =  29.237  Annual  Evap(cm)  =  30.964 


YEAR  =  600 

CAC03(gCa/cm2)  =  0.13084E-04  0.39974E-04  0.59466E-03  0.17612E-01  0.13134E-01 
CAC03(gCa/cm2)  =  0.78798E-02  0.39994E-02  0.12391E-02  0.40816E-03  0.18832E-03 

CAS04 (gS/cm2)  =  O.OOOOOE+00  O.-OOOOOE+OO  O.OOOOOE+00  O.OOOOOE+00  O.OOOOOE+00 
CAS04 (gS/cm2)  =  O.OOOOOE+00  0.19832E-03  0.93836E-03  0.19800E-02  0.45282E-03 

Bulk  Density (g/cm3)  =  1.4400  1.4400  1.4401  1.4444  1.4433 

Bulk  Density (g/cm3)  =  1.4420  1.4411  1.4408  1.4412  1.4403 

Leach (cm)  =  0.6681E+01  Leacca (gCa/cm2)  =0.47531E-03  Leachs (gS/cm2)  =0.17596E-03 
Annual  Rain (cm)  =  21.983  Annual  Evap(cm)  =  21.989 


YEAR  =  800 

CAC03(gCa/cm2)  =  0.40850E-05  0.56173E-04  0.20961E-03  0.24291E-01  0.17550E-01 
CAC03(gCa/cm2)  =  0.10510E-01  0.53865E-02  0.16162E-02  0.49783E-03  0.21882E-03 

CAS04 (gS/cm2)  =  O.OOOOOE+00  O.OOOOOE+00  O.OOOOOE+00  O.OOOOOE+00  O.OOOOOE+00 
CAS04 (gS/cm2)  =  O.OOOOOE+00  0.51036E-03  0.85479E-03  0.29788E-02  0.79254E-03 

Bulk  Density (g/ cm3)  =  1.4400  1.4400  1.4401  1.4461  1.4444 

Bulk  Density (g/ cm3)  =  1.4426  1.4416  1.4409  1.4417  1.4405 

Leach (cm)  =  0.6681E+01  Leacca (gCa/cm2)  =0.47531E-03  Leachs (gS/cm2)  =0.17596E-03 
Annual  Rain (cm)  =  30.808  Annual  Evap(cm)  =  30.486 


YEAR  =  1000 

CAC03(gCa/cm2)  =  O.OOOOOE+00  0.69365E-04  0.30493E-03  0.30038E-01  0.22281E-01 
CAC03(gCa/cm2)  =  0.13423E-01  0.65728E-02  0.20215E-02  0.59955E-03  0.22829E-03 

CAS04 (gS/cm2)  =  O.OOOOOE+00  O.OOOOOE+00  O.OOOOOE+00  O.OOOOOE+00  O.OOOOOE+00 
CAS04 (gS/cm2)  =  O.OOOOOE+00  0.20971E-03  0.12965E-02  0.42528E-02  0.92806E-03 

Bulk  Density (g/cm3)  =  1.4400  1.4400  1.4401  1.4475  1.4456 

Bulk  Density (g/ cm3)  =  1.4434  1.4418  1.4412  1.4424  1.4406 

Leach(cm)  =  0.6681E+01  Leacca (gCa/cm2)  =0.47531E-03  Leachs (gS/cm2)  =0.17596E-03 
Annual  Rain (cm)  =  31.399  Annual  Evap(cm)  =  29.940 


FINAL  SOIL  PROFILE  CONCENTRATIONS 
See  last  year  above  for  other  final  soil  properties 
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Table  2  (Cont'd). 


4/18/94  10:37 


Macintosh  HD:MPW;CaData 


Hori 

Mol SCO (g/cm2) 

H  activity 

Ca (gCa/cm2) 

1 

0.81111E+00 

0.71021E-08 

0.15696E-04 

2 

0.89038E+00 

0.67547E-08 

0.43346E-04 

3 

0.14261E+01 

0.89597E-08 

0.91333E-04 

4 

0.67716E+00 

0.12201E-07 

0.69966E-04 

5 

0.56160E+00 

0.15032E-07 

0.76210E-04 

6 

0.56160E+00 

0.23877E-07 

0,23684E-03 

7 

0.56160E+00 

0.28090E-07 

0,32912E-03 

8 

0.56160E+00 

0.28948E-07 

0,32946E-03 

9 

0.56160E+00 

0.29846E-07 

0.32981E-03 

10 

0.56160E+00 

0.30658E-07 

0.33013E-03 

S04 (gS/cm2) 
0.14565E-05 
0.64568E-05 
0.24335E-04 
0.33322E-04 
0.41713E-04 
0.17419E-03 
0.24847E-03 
0.24831E-03 
0.24814E-03 
0.24798E-03 


Torain(cm)  =  0.286059E+05  Toevap(cm)  =  0.286064E+05 

Residual  Dust  Ca(gCa/cm2)  =  0.2360E-05  Residual  Dust  S04(gS/cm2)  =  0.2190E-06 


data  necessary  to  reinitialize  the  model  for  runs 
with  new  climates  or  other  system  drivers.  One 
could,  for  example,  simulate  alternating  Pleistocene 
and  Holocene  climates. 

Rainmodule  output 

For  convenience,  the  current  version  of 
Rainmodule  prints  to  the  computer  screen.  One 
can  get  a  hard  copy  by  selecting  "Print  Window" 
to  sent  the  output  to  the  on-line  printer.  Alterna¬ 
tively,  the  user  could,  if  desired,  change  the  pro¬ 
gram  to  "Write"  directly  to  a  file  as  was  done  in 
CALGYP. 

This  program  prints  total  rainfall  and  total  rain 
events  for  a  fixed  simulation  of  1000  years  in  100- 
year  blocks  (Table  3).  The  average  annual  rainfall 
should  be  close  to  the  long-term  mean  for  the  site. 
In  this  particular  case  for  Tucson,  calculated  aver¬ 
age  annual  rainfall  was  28.6  cm  compared  to  28.4 
cm  for  the  long-term  mean  (Table  3)  .  Note  that  the 
total  rainfall  in  1000  years  (28,605.9  cm)  is  identical 
for  both  CALGYP  (T  able  2)  and  Rainmodule  (Table 
3)  because  the  same  Seed  (123)  was  used  for  the 
random-number  generators  in  both  programs. 
Using  a  different  Seed  of  1  yields  44,833  rain 
events  for  a  total  of  28,675.2  cm  of  ram  in  1000 
years,  which  is  similar  but  not  identical  to  the 
previous  run  with  a  Seed  of  123  (Table  3). 


PROGRAM  VALIDATION 

CALGYP  maintains  massbalances  with  respect 
to  water,  calcium,  and  sulfate.  Users  should  verify 
for  themselves  that  the  model  accurately  main¬ 
tains  these  balances.  CALGYP  prints  the  necessary 
output  to  make  these  calculations. 


Water  balance  is  given  by 

(Soil  Water)jx\itial  +  Rain  =  (Soil  Water)fjnal  + 
Evapotranspiration  +  Leachate. 

For  the  1000-year  simulation  (Table  2),  these  quan¬ 
tities  (cm  of  water)  are 

14.4  +  28,605.9  =  7.2  +  28,606.4  +  6.7 
28,620.3  =  28,620.3. 

For  calcium,  the  major  sink  in  the  soil  profile  is 
CaCOa  which  has  accumulated  7.553  x  10~2  g  Ca 
cm“2  over  the  1000-year  simulation  (T able  2) .  Other 
sinks  include  CaS04  (8.36  x  10“^  g  Ca  cm-2),  A 
soluble  Ca*  (+1.77  x  10~3  g  Ca  cm-2),  leachate  Ca 
(4.8  X 10“^  g  Ca  cm-2),  and  residual  dust  Ca  (2x10“ 
^  g  Ca  cm“2).  The  soil  profile  increased  in  Ca 
content  over  the  1000-year  simulation  by  8.614  x 
10-2  g  (3a  cm-2,  which  is  equal  to  the  annual 
atmospheric  input  of  8.614  x  10“^  g  Ca  cm-2 
year-1  (Table  2)  for  1000  years. 

For  sulfate,  the  major  sink  is  CaS04*2H20  which 
has  accumulated  6.687  x  10-3  g  g  cm-2  over  the 
1000-year  simulation  (Table  2).  Other  sinks  in¬ 
clude  A  soluble  SO4  (+1.131  x  10-3  g  5  cm-2), 
leachate  SO4  (1.76  x  10-^  g  S  cm-2),  and  residual 
dust  SO4  (2  X  10-2  g  S  cm-2).  The  soil  profile 
increased  in  SO4  by  7.994  x  10-3  g  g  cm-2  which  is 
equal,  within  roimdoff  error,  to  the  annual  SO4 
dust  input  rate  of  7.993  x  10-^  g  S  cm-2  year-i 
(Table  2)  for  1000  years. 

For  this  particular  simulation,  only  dust  chemi- 


*The  change  (A)  in  soluble  Ca  between  the  final  (1000 
yr)  and  initial  (0  yr)  of  the  simulation. 
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Table  3.  Sample  input  and  output  of  the  rain 
module  program  for  a  1000-year  simulation  for 
the  Tucson  site. 


Select  Site: 

1  =  El  Paso,  TX  (21.6an) 

2  =  Albuq.,  NM  (21.1cm) 

3  =  Clayton,  NM  (37.8cm) 

4  =  Roswell,  NM  (31.6cm) 

5  =  Yuma,  AZ  (8.5cm) 

6  =  Phoenix,  AZ  (18.9cm) 

7  =  Tucson,  AZ  (28.4cm) 

7 

Title? 

Tucson  Simulation 

Fractional  Change  in  Rain  for  Drier  Climate. 
Default  =  0.00,  no  change. 


0.0 

Seed  for  Random  Number  Generator? 


123 

100  yrs 
100  yrs 
100  yrs 
100  yrs 
100  yrs 
100  yrs 
100  yrs 
100  yrs 
100  yrs 
100  yrs 
1000  yrs 


chemical  concentrations  are  also  specified,  then 


Rain (cm)  Rain  Events 

3049.7 

4554 

2871.1 

4537 

2892.7 

4442 

2885.6 

4335 

2879.0 

4479 

2788.4 

4335 

2840.2 

4386 

2712.4 

4273 

2867.2 

4393 

2819.7 

4395 

28605.9 

44129 

■re  specified 

(Table  2).  When  rain 

Ca  (or  SO4-S)  (g  cm-2)  =  Ca  (or  SO4  -  S) 
(mg  L-i)  X  Rain  (cm)  x  1.0  x  10^ 


must  be  added  to  dust  inputs  to  determine  total 
system  chemical  input. 

The  stochastic  rain  model  used  in  CALGYP 
accurately  predicts  the  mean  annual  rainfall  and 
the  variability  in  this  quantity  for  the  seven  South¬ 
western  desert  sites  (Marion  et  al.  1985).  The 
CALGYP  model  using  current  climate  typically 
predicts  a  shallower  depth  of  CaCOs  deposition 
than  is  found  in  most  Southwestern  desert  soils. 
This  is  not  surprising  because  most  CaCOs  is 
believed  to  have  formed  xmder  earlier,  "wetter" 
Pleistocene  climates  (Gile  et  al.  1981,  Marion  et  al. 
1985,  McFadden  and  Tinsley  1985).  CALGYP  is 
compatible  with  field  soils  if  one  assumes  that 
most  pedogenic  CaCOa  formed  during  a  cool-wet 
Pleistocene  climate  (Marion  et  al.  1985).  CALGYP 
predicts  that  the  time  required  for  complete  plug¬ 
ging  of  soil  profiles  with  CaCOa  requires  » 10,000 
years  in  agreement  with  independent  evidence 
(Gile  et  al.  1981;  Marion  et  al.  1985).  CALGYP 
correctly  predicts  the  deeper  deposition  of  the 


more  soluble  mineral,  gypsum,  in  the  soil  profile 
relative  to  calcite  (Fig.  8).  The  dependence  of  soil 
pFI  on  sulfate  concentration  (eq  18)  was  validated 
recently  with  field  data  from  a  cold  dry  site  in 
Alaska  and  a  hot  dry  Mojave  Desert  site  (Marion 
and  Schlesinger  in  press).  Other  aspects  of  CALGYP 
have  been  validated  previously  (Marion  et  al.  1985, 
Marion  and  Schlesinger  in  press). 


ALTERED  CLIMATE  AND  CHEMISTRY 
Climate 

CALGYP  is  structured  to  alter  both  tempera¬ 
ture  and  rainfall  for  the  seven  Southwestern  sites 
currently  in  the  model.  There  are,  however,  limita¬ 
tions  on  these  alterations. 

Temperature  largely  controls  output  of  water 
through  the  evapotranspiration  process  (Fig.  4 
and  6).  To  keep  temperature  within  the  range  of 
data  used  to  parameterize  the  model,  mean 
monthly  temperatures  should  be  between  5  and  32 
°C  (Fig.  4)  and  mean  annual  temperatures  should 
be  between  12  and  23  °C  (Fig.  6).  Any  temperature 
alteration  should  ideally  be  limited  to  these  ranges. 
However,  even  in  the  current  climate  model,  tem¬ 
peratures  occasionally  faU below  5  °C  in  the  winter 
months  for  Albuquerque,  Clayton,  and  Roswell 
(App.  A),  but  mean  monthly  temperatures  never 
fall  below  0  °C.  In  order  to  keep  temperatures 
within  the  ideal  range,  one  could,  for  example, 
either  decrease  mean  monthly  temperature  by  5 
°C  for  Yuma,  Phoenix,  or  Tucson  or  raise  mean 
monthly  temperature  by  5  °C  for  El  Paso,  Albu- 


Figure  8.  Accumulation  of  calcite  and  gypsum  in  a 
1000-year  simulation  for  the  Tucson  site  (data  in 
Table!). 
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querque,  Clayton,  or  Roswell,  but  not  vice-versa 
(App.  A).  The  temperature  change  is  specified  in  a 
screen  prompt  under  "Altered  Climate." 

Rain  patterns  can  be  altered  by  simply  chang¬ 
ing  the  data  in  Detrain  for  the  deterministic  model 
or  by  changing  the  frequency  distributions  for  the 
stochastic  model.  To  alter  the  rain  pattern  in  the 
stochastic  model  requires  explicitly  changing  "sea¬ 
sonal"  patterns.  In  the  current  model,  the 
interarrival  days  (variable  "Inter,"  App.  A)  for 
Tucson  in  the  winter  is  given  by 

0, 1, 2, 3, 4,  7, 10, 15, 25, 40, 57. 

Changing  these  days  to: 

0, 1, 2, 3, 4, 5, 6,  8, 10, 12, 15 

while  holding  the  corresponding  probabilities  in 
"Freq"  constant  has  the  result  of  decreasing  the 
interarrival  time  between  rain  events,  and  thereby 
increasing  the  number  of  winter  rainfall  events. 
Similarly,  changing  the  rainfall  amoimts  (variable 
"R")  for  Tucson-winter  from 

0.00, 0.25, 0.50, 0.75, 1 .00, 1.25, 1 .50, 1 .75, 2.00, 
3.00, 7.52 

to 

0.00, 0.25, 0.50, 0.75, 1 .00, 1 .25, 1 .50, 2.00, 3.00, 
5.00, 9.50 

while  holding  the  corresponding  probabilities  in 
"D"  constant  results  in  increasing  the  average 
rainfall  amount/event.  Together,  these  two 
changes  in  interarrival  time  and  rainfall  amoimts 
increase  the  mean  annual  rainfall  for  Tucson  from 
28.6  cm  (Tables  2  and  3)  to  43.4  cm,  a  51 .6  %  change 
with  all  of  the  increase  occurring  in  the  winter 
months.  Most  of  this  change  is  caused  by  increas¬ 
ing  the  number  of  rainfall  events  from  44.1  year-i 
(Table  3)  to  60.4  yeam^  (37.0%  increase)  and  the 
remainder  (14.6%)  is  caused  by  increased  rain 
amounts/event.  By  trial-and-error,  one  can  ma¬ 
nipulate  the  stochastic  model  probability  distribu¬ 
tions  to  increase  or  decrease  rain  dates,  rainfall 
amounts,  or  alter  seasonal  distributions. 
Rainmodule  is  included  with  CALGYP  to  facili¬ 
tate  these  trial-and-error  calculations. 

Provided  temperatures  remain  within  the  range 
of  model  parameterization  (see  previous  discus¬ 
sion),  the  user  can  easily  add  new  Southwestern 
desert  sites  to  CALGYP.  Adding  new  sites  to 


CALGYP  from  outside  the  desert  Southwest  will 
require  developing  a  new  evapotranspiration  al¬ 
gorithm.  One  needs  to  relate  actual  evapotranspi¬ 
ration  to  the  soil-vegetative  system  of  interest 
(e.g.,  deserts,  forests,  grasslands)  (Fig.  5).  Such  an 
exercise  is  not  trivial,  but  neither  is  it  a  major 
problem.  In  most  cases,  the  type  of  information 
needed  to  develop  these  relationships  is  available 
in  the  literature.  To  extrapolate  this  model  to  sites 
where  seasonal  freezing  is  important  will  require 
changing  how  seasons  are  dealt  with.  For  ex¬ 
ample,  winter  might  accumulate  water  as  snow, 
allowing  for  a  single  leaching  event  with  some 
fraction  of  winter  snow  in  the  spring.  Because 
CALGYP  is  structured  by  seasons,  incorporating 
frozen  seasons  should  not  be  a  major  problem. 

Chemistry 

InEquil,  chemical  equilibrium  is  maintained  for 
the  CaS04°  ion-pair  and  the  precipitation-disso¬ 
lution  of  calcite  and  gypsum.  CALGYP  uses  a 
sequential  approach  to  solving  the  chemical  equi¬ 
librium  equations.  This  approach  is  both  simple 
and  flexible.  FREZCHEM,  a  chemical  thermody¬ 
namic  model  for  aqueous  solutions  at  subzero 
temperatures,  uses  ^is  approach  and  deals  with 
ice  formation  and  the  precipitation-dissolution  of 
15  minerals  (Marion  and  Grant  1994).  To  add  new 
reactions  means  adding  new  chemical  algorithms 
to  the  present  sequence  (See  Equil).  Depending  on 
what  the  new  reactions  involve,  the  convergence 
criteria  may  need  to  be  changed.  At  present  the 
model  iterates  until  the  Ca2+  ion  concentrations  in 
successive  iterations  agree  to  within  ±  1%.  Other 
solution  species  or  a  suite  of  species  could  be  used 
to  test  for  mathematical  convergence. 

The  upper  limit  of  solution  concentrations  that 
CALGYP  can  handle  is  set  by  the  range  of  validity 
of  the  Davies  equation  (eq  14),  which  is  approxi¬ 
mately  0.1  M  (Davies  1962).  Because  of  this  limita¬ 
tion,  adding  a  very  soluble  salt  such  as  NaCl  is  not 
currently  feasible.  The  solubility  of  NaCl  at  25  °C 
is  approximately  6.1  M  (Marion  and  Grant  1994). 
To  develop  a  model  capable  of  working  athigh  salt 
concentrations  can  be  done  using  the  Pitzer  Equa¬ 
tions  (Plummer  et  al.  1988,  Spencer  et  al.  1990, 
Pitzer  1991,  Marion  and  Grant  1994). 
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PROGRAM  CALGYP 

C  This  is  CALGYP.  The  primary  function  is  to  calculate  CaC03  and 
C  CaS04  deposition  in  soils.  This  program  will  run  using  current 
C  climatic  drivers  for  El  Paso,  Albuquerque,  Clayton,  Roswell,  Yuma, 

C  Phoenix,  and  Tucson.  The  model  will  also  simulate  altered  climates. 

C  You  may  specify  a  drier  climate  or  a  change  in  temperature  directly; 

C  see  following  input  statements.  See  technical  report  for 
C  restrictions  on  temperature  change.  To  specify  a  wetter  climate 
C  requires  manipulating  the  cumulative  probability  distributions  for 
C  rain  dates  or  rainfall  amounts  for  the  appropriate  season  and  site. 

C  See  accompanying  technical  report  for  suggestions  on  how  this  might 
C  be  done. 

REAL  Midmc,  Me,  Ll,Isum,  Leach,  Leacca,  Leachs 
REAL  Mpotev (12) ,Midmoi (10) , Moist (10) , Lowerl (10) 

REAL  ICO2(10) ,IBd(10) ,IHs(10)  ,Loss(10) 

REAL* 4  Rnum, rand 

INTEGER  Site, Climat, Hori, Year, Years, Pint, S, Sumfre, Days (12) , Seed 
DIMENSION  Thic)c  (10)  ,  Bd(lO),  BarOl(lO),  Barl5(10), 

X  Ca(lO),  C02Mul(12),  S04 (10) ,  Temp (7, 12), 

X  Upperl(lO),  Whc(lO) ,Tis(12) ,  Dpotev(12), 

X  Availw(lO),  Poresp(lO) 

CHARACTER*50  Title 

DOUBLE  PRECISION  S04s, CaC03, CaS04, Cas, Ka, Hs 
DOUBLE  PRECISION  Moisco, Horizo, Gain 

COMMON  /  Precipitation  /  Rain (200),  Rainda (200) , Simfre, Site, 

X  Delrai, Numb, Rnum, Inter (7, 3, 11) , Freq (7, 3,11),R(7,3,11), 

X  D(7,3,ll) 

COMMON  /  Chemistry  /  S04s (10) , CaC03 (10) , CaS04 (10) , Cas (10) ,Dhc, 

X  C02 (10) ,Hs (10) ,Ka(6) , Moisco (10) , Horizo (10) , Gain (10) 

OPEN  (UNIT  =  2,  FILE  =  "CaData") 

C  Read  in  the  program  data 

PRINT*,  'Select  Site:' 

PRINT*,  '  1  =  El  Paso,  TX  (21.6cm)' 

PRINT*,  '  2  =  Albuq.,  NM  (21.1cm)' 

PRINT*,  '  3  =  Clayton,  NM  (37.8cm)' 

PRINT*,  '  4  =  Roswell,  NM  (31.6cm)' 

PRINT*,  '  5  =  Yuma,  AZ  (8.5cm)' 

PRINT*,  '  6  =  Phoenix,  AZ  (18.9cm)' 

PRINT*,  '  7  =  Tucson,  AZ  (28.4cm)' 

READ*,  Site 

Deitem  =  0.0 
Delrai  =  0.0 

PRINT*,  'Climate  Option?  Current  Climate  =  1,  Altered  Climate  =  2' 
READ*,  Climat 
IF  (Climat  .eq.  2)  THEN 

PRINT*, 'Enter  Delta  Temperature=  (°C, default=0 . 0, no  change)' 
READ*,  Deitem 

PRINT*,  'Enter  Fractional  Change  in  Rainfall  Amount  during' 
PRINT*,  '  Drier  Climate  =  (default  =0.0,  no  change)' 

READ*,  Delrai 

END  IF 

C  Days/Month 

DATA  Days  /31,  28,  31, 30,  31, 30,  31,  31, 30, 31,  30,  31/ 
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C  Monthly  Temperatures 

•DATA(  (Tempd,  J)  ,  J=l,12)  ,1=1,7)  / 

X  7.1,9.6,13.1,17.6,22.3,27.1,27.8,26.8,23.8,18.1,11.3,7.3,  El  Paso 

X  1.4,4.3,8.1,12.7,17.7,23.1,25.2,24.1,20.3,13.8,6.6,1.8,  Albuq 

X  0.6,2.4,5.2,10.7,15.6,20.9,23.3,22.4,18.4,12.8,5.6,2.1,  Clayton 

X  4.2,6.6,10.6,15.3,20.1,24.8,26.2,25.4,21.6,15.4,8.8,4.5,  Roswell 

X  12. 8, 15. 1,17. 8, 21. 4, 25. 2,  29. 8,  33. 4,  32. 9, 29. 8, 23. 5, 17. 2, 13. 3, Yuma 

X  11. 1,13. 2, 15. 9, 20. 0,24. 6, 29. 7,  32. 8,  31. 8, 28. 8, 22. 2, 15. 6, 11. 6, Phoenix 

X  10.2,11.8,14.3,18.1,22.6,27.8,30.0,28.9,26.7,20.8,14.6,10.8/Tucson 

C  Read  in  Stochastic  Rain  Model  Parameters 

DATA  (((D(I,J,K),K=1,11),J=1,2),I=1,7)  / 

X  .000, .571, .745, .827, .837, .888, .918, .939, .959, .990,1.000, 

X  .000, .526, .691, .773, .835, .866, .897, .938, .948, .979,1.000, 

X  .000, .635, .832, .920, .964, .971, .986, .993,1.000,1.05,1.10, 

X  .000,  .577,  .740,  .827, .904, .952, .962, .981, .991,1.00,1.05, 

X  .000, .586, .811, .928, .955, .982, .991,1.00,1.05,1.10,1.15, 

X  .000,  .371,  .596,  .728, .794, .860, .880, .933, .966, .986,1.000, 

X  .000, .570, .756, .884, .907, .942, .965,1.00,1.05,1.10,1.15, 

X  .000, .405, .628, .719, .769, .827, .868, .901, .975, .983,1.000, 

X  .000, .512, .780, .865, .889, .938, .962, .986,1.00,1.05,1.10, 

X  .000, .625, .875,1.000,1.05,1.10,1.15,1.20,1.25,1.30,1.35, 

X  .000, .394, .567, .711, .817, .875, .923, .952, .971, .990,1.000, 

X  .000, .600, .800, .933,1.000,1.1,1.2,1.3,1.4,1.5,1.6, 

X  .000,  .473,  .661, .732, .803, .830, .875, .902, .947, .  983,1.000, 

X  .000, .550, .850, .950,1.000,2.0,3.0,4.0,5.0,6.0,7.0  / 

DATA  ( (D (I, 3,K) ,K=1,11) ,1=6,  7)  / 

X  .000, .642, .755, .793, .868, .925, .963, .982,1.000,1.5,2.0,  PH, SUM 

X  .000, .426, .539, .617, .730, .765, .843, .878, .913, .991,1.000  /  TU,SUM 

DATA  (((R(I, J,K),K=1,11),J=1,2),I=1,7)  / 

X  0.00, .25, .50, .75,1.00,1.25,1.50,1.75,2.00,3.00,3.07, 

X  0.00, .25, .50, .75,1.00,1.25,1.50,1.75,2.00,4.00,5.59, 

X  0.00, .25, .50, .75,1.00,1.25,1.50,2.00,2.64,3.00,5.00, 

X  0.00, .25, .50, .75,1.00,1.25,1.50,2.00,3.00,4.45,10.0, 

X  0.00, .25, .50, .75,1.00,1.25,2.00,4.83,5.00,10.0,15.0, 

X  0.00,  .25,  .50,  .75,1.00,1.50,2.00,3.00,4.00,5.00,5.36, 

X  0.00, .25, .50, .75,1.00,1.25,1.50,1.75,5.00,10.0,15.0, 

X  0.00, .25, .50, .75,1.00,1.50,2.00,3.00,5.00,7.00,10.72, 

X  0.00, .25, .50, .75,1.00,1.25,1.50,3.00,4.45,10.0,16.0, 

X  0.00,  .25,  .50, .58,1.00,1.25,1.50,1.75,2.00,3.00,3.07, 

X  0.00, .25, .50, .75,1.00,1.25,1.50,1.75,2.00,3.00,5.03, 

X  0.00,  .25, .50, .75, .97,1.25,1.50,1.75,2.00,3.00,3.07, 

X  0.00, .25, .50, .75,1.0,1.25,1.5,1.75,2.0,3.0,7.52, 

X  0.00, .25, .50, .75,1.02,2.0,3.0,4.0,5.0,6.0,7.0  / 


DATA  ( (R(I,3,K) ,K=1,11) ,1=6,7)  / 

X  0.00, .25, .50, .75,1.00,1.50,2.00,3.50,4.22,5.00,10.0,  PH, SUM 

X  0.00, .25, .50, .75,1.00,1.25,1.50,1.75,2.00,3.00,3.23  /  TU, SUM 

DATA  ( ( (Freq(I, J,K) ,K=1,11) , J=l,2) ,1=1,7)  / 

X  0.00, .388, .449, .510, .531, .643, .765, .796, .888, .969,1.000,  EP,WINT 

X  0.00, .402, .454, .567, .680, .750, .853, .902, .923, .974,1.000,  EP,SUM 

X  0.00, .328, .447, .522, .589, .664, .768, .843, .910, .977,1.000,  AL,WINT 

X  0.00, .396, .481, .594, .679, .783, .840, .963, .982, .991,1.000,  AL, SUM 

X  0.00, .291, .409, .427, .472, .617, .762, .844, .944, .989,1.000,  CL,WINT 

X  0.00, .404, .530, .623, . 716, . 848, . 901, .961, .987, .994,1.000,  CL, SUM 

X  0.00, .430, .488, .569, .616, .744, .825, .918, .988,1.000,1.05,  RS,WINT 

X  0.00, .380, .430, .521, .554, .703, .802, .918, .984, .992,1.000,  RS,SUM 


EP,WINT 

EP,SUM 

AL,WINT 

AL,SUM 

CL,WINT 

CL, SUM 

RS,WINT 

RS,SUM 

YU,WINT 

YU, SUM 

PH,WINT 

PH,SPR 

TU,WINT 

TU, SPR 


EP,WINT 

EP,SUM 

AL,WINT 

AL, SUM 

CL,WINT 

CL, SUM 

RS,WINT 

RS,SUM 

YU,WINT 

YU, SUM 

PH,WINT 

PH, SPR 

TU,WINT 

TU, SPR 
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X  0.00,  .316,  .417, .430, .4  93, .582, .  671,  . 78  5, . 950,  . 988, 1 . 000,  YU,WINT 

X  0.00, .111, .333, .777,1.000,1.05,1.1,1.2,1.3,1.4,1.5,  YU, SUM 

X  0.00, .476, .515, .583, .641, .719, .787, .865, .923, .972,1.000,  PH,WINT 

X  0.00, .250, .313, .438, .501, .564, .689, .814, .939,1.00,1.05,  PH,SPR 

X  0.00, .447, .552, .596, . 631, . 745, . 815, . 885, .955, .973,1.000,  TU,WINT 

X  0.00, .222, .278, .334, .390, .501, .612, .779, .890, .946,1.000  /  TU,SPR 

DATA  ((Freq{I,3,K),K=l,ll),I=6,7)  / 

X  0.00, .212, .347, .482, .559, .713, .790, .886, .963,1.00,1.05,  PH, SUM 

X  0.00, .491, .640, .728, .781, .886, .965, .991,1.00,1.1,1.2  /  TU,SUM 

DATA  (((Inter(I,J,K),K=l,ll),J=l,2),I=l,7)  / 

X  0,1,2,3,4,7,10,15,20,30,57, 

X  0,1,2,3,4,7,10,15,20,30,66, 

X  0,1,2,3,4,7,10,15,20,30,40, 

X  0,1,2,3,4,7,10,15,20,30,32, 

X  0,1,2,3,4,7,10,15,25,40,50, 

X  0,1,2,3,4,7,10,13,15,20,21, 

X  0,1,2,3,7,10,15,25,40,64,100, 

X  0,1,2,3,4,7,10,15,25,40,48, 

X  0,1,2,3,4,7,10,15,40,70,111, 

X  0,20,40,60,121,150,160,170,180,190,200, 

X  0,1,2,3,4,7,10,15,25,40,52, 

X  0,1,5,10,15,25,40,50,60,77,100, 

X  0,1,2,3,4,7,10,15,25,40,57, 

X  0,1,3,5,7,10,25,30,40,60,78  / 

DATA  ((Inter(I,3,K),K=l,ll),I=6,7)  / 

X  0,1,2,3,4,7,10,15,20,24,50,  PH, SUM 

X  0,1,2,3,4,7,10,15,17,40,50  /  TU,SUM 

C  Basic  Soil  &  Atm  Parameters 

C  Horizon  Thickness  (cm)  and  Bulk  Density  (g/cc) 

DATA  Thick  /  10*10.0  / 

DATA  Bd  /  10*1.44  / 

C  Horizon  Soil  Water  Cone.  (%)  at  0.1  Bar,  15  Bar,  and  Initial 
DATA  BarOl  /  10*. 122  / 

DATA  Barl5  /  10*. 039  / 

DATA  Moist  /  10*0.10  / 

C  Horizon  CaC03  (%Ca)  and  Soluble  Ca  (%) 

DATA  CaC03  /  10*0.0  / 

DATA  Ca  /  10*6e-7  / 

C  Horizon  H  Activity  and  C02  Cone.  (ATM) 

DATA  Hs  /  10*le-8  / 

DATA  C02  /3.79e-4,6.87e-4,9.76e-4,1.28e-3,1.60e-3, 

X  1. 91e-3,2.14e-3,2.27e-3,2.41e-3,2. 54e-3/ 

C  Monthly  C02  Multiplier 

DATA  C02Mul  /  2*1.00,8*1.588,2*1.00  / 

C  Horizon  CaS04  (%S)  and  Soluble  S04  (%S) 

DATA  CaS04  /  10*0.0  / 

DATA  S04  /10*1.0e-6/ 

C  Dust-Ca  (g/cm2/day) ,  Rain-Ca  (mg/1),  Rain-S04-S  (mg/1), 

C  Dust-S04-S  (g/cm2/day) 

Dustca  =  2.36e-7 
Precca  =  0.00 
Precs  =  0.00 
Dusts  =  2.19e-8 

C  Read  In  Basic  Model  Run  Parameters 

PRINT*,  'Title?' 

READ (*,10)  Title 
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10  FORMAT (A50) 

PRINT*,  'Nvunber  of  Soil  Horizons  (Max=10)?' 

READ*,  Hori 

PRINT*,  'Years  to  Simulate?' 

READ*,  Years 

PRINT*,  'Print  Interval?' 

READ*,  Pint 

PRINT*,  'Deterministic  (1)  or  Stochastic  (2)  Rainfall  Model?' 
READ*,  S 

IF  (S  .eq.  1)  GOTO  15 

PRINT*,  'Seed  for  Random  Number  Generator?' 

READ*,  Seed 
Rnum  =  rand (Seed) 

C  Convert  Input  to  Program  Units  {%  to  g/cm2) 

15  DO  20  I  =  1,  Hori 

Moisco (I) =Thick (I) *Bd(I) *Moist (I) 

Lowerl  (I)  =Thic)c  (I)  *Bd(I)  *Barl5  (I) 

Upperl (I)=Thick(I)*Bd(I)*Bar01 (I) 

Whc (I) =Upperl  (I) -Lowerl (I) 

Diff=Whc(I)*.546 

Midmoi (I) =Lowerl (I) +Diff 

CaC03 (I) =CaC03 (I) *Thick (I) *Bd(I) 

Ca(I)=Ca(I) *Thick(I) *Bd(I) 

CaS04 (I)=CaS04 (I) *Thick (I) *Bd{I) 

S04 (I)=S04 (I)*Thick(I)*Bd(I) 

20  CONTINUE 

C  Calculate  Monthly  Potential  Evaporation  using 
C  Thornthwaite's  Equation  and  Adjust  to  Pan  Evaporation 

I sum  =  0.0 
Totalt  =  0.0 
DO  30  J  =  1,  12 

Tis(J)  =  (TempCSite, J)/5.0)**(1.514) 

Isum  =  Isum  +  Tis(J) 

Totalt  =  Totalt  +  Temp (Site,  J) 

30  CONTINUE 

Avert  =  Totalt/12.0 

Rat  =  (0.1278*EXP (.2802* (Avert+Deltem) )+225. 9) / 

X  (0.1278*EXP (.2802*Avert)+225.9) 

A=(6.75e-7* (lsum**3.0) ) - (7 . 71e-5* (Isum**2 . 0) )+ 

X  (l,792e-2*Isum+. 49239) 

DO  40  J  =  1,  12 

Mpotev(J)  =  1 . 6* ( (10 . 0*Temp (Site, J) /Isum) **A) *Days ( J) /30 . 0 
Dpotev(J)  =  Mpotev(J)  /  Days(J) 

IF  (J  .le.  7)  THEN 

Y  =  17.07*EXP(-.1309*Temp(Site, J) )+1.91 

ELSE 

Y  =  13. 67*EXP(-.1300*Temp (Site, J) ) +1.35 

END  IF 

Dpotev(J)  =  Y*Dpotev ( J) *Rat 
Mpotev(J)  =  Y*Mpotev ( J) *Rat 
Temp (Site, J)  =  Temp (Site, J) +Deltem 
40  CONTINUE 

C  Print  the  Initial  State  of  the  System 

WRITE  (2,42)  Title 
42  FORMAT  (A50) 

WRITE  (2,*)  'SITE  =',  Site 

IF  (S  .eq.  1)  WRITE  (2,*)  'Deterministic  Rainfall  Model' 
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IF  (S  .eq.  2)  WRITE  (2,*)  'Stochastic  Rainfall  Model  with  Seed  =', 
X  Seed 

WRITE  (2,*) 

WRITE  (2,*)  ’  CLIMATIC  SUMMARY’ 

WRITE  (2,*)  '  Month  Temp(C)  Pot.Evap. (g/cm2) ' 

DO  50  I  =  1,  12 

WRITE  (2,45)  I,  Temp (Site, I),  Mpotev(I) 

45  FORMAT(10X,I2,6X,F4.1,8X,F5.2) 

50  CONTINUE 

WRITE  (2,*) 

WRITE  (2,*)  'INITIAL  SOIL  PROFILE' 


WRITE  (2,*) 

'Hori  Thick 

BD  C02 

H  CaC03 

:  Ca 

CaS04 

S04' 

WRITE  (2,*) 

'  cm 

g/cm3  atm 

gCa/cm2 

:  gCa/cm2 

gS/cm2 

gS/cm2 ' 

DO  60  I  =  1, 

Hori 

WRITE  (2,55)  I,  Thic)c{I),  Bd(I),  C02(I),  Hs(I),  CaC03(I), 
X  Ca(I),  CaS04(I),  S04(I) 

55  FORMAT (IX, I2,F6.1,F6.2, 6E11,3) 

60  CONTINUE 

WRITE  (2,*) 

WRITE  (2,*)  'ATMOSPHERIC  CHEMICAL  CONDITIONS' 

Dca=Dustca*365 . 0 
Ds=Dusts*365 . 0 

WRITE  (2,*)  ' Dustca (gCa/cm2/yr)  Precca (mgCa/1) 

X  Dusts (gS/cm2/yr)  Frees (mgS/1) ' 

WRITE  (2,65)  Dca,  Precca,  Ds,  Frees 
65  FORMAT  (2X,E12.4,F12.2,11X,E12.4,F12.2) 

C  Initialize  the  Program  Parameters 

DO  70  I  =  1,  Hori 

Ca(I)  =  Ca(I)/40.08 
S04  (I)  =  S04 (I) /32.064 
IC02 (I)  =  C02 (I) 

IBd(I)  =  Bd(I) 

IHs(I)  =  Hs(I) 

70  CONTINUE 

Year  =  1 
Leach  =  0 
Leacca  =  0 
Acumca  =  0 
Leachs  =  0 
Accums  =  0 
Torain  =  0 
Toevap  =  0 
Inte  =  0 
Cement  =  0 
100  Rdays  =  0 

Tdays  =  Days(l) 

J=1 

1=1 

CALL  CONSTANTS (Ka,  Temp (Site, 1) , Dhc) 

Tevap  =  0 
Train  =  0 
Day  =  0 
Midmc  =  0 
Me  =  0 
LI  =  0 
U1  =  0 

DO  110  K  =  1,  Hori 

LI  =  Ll  +  Lowerl (K) 

U1  =  U1  +  Upperl (K) 
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Mldmc  =  Midmc  +  Midmoi (K) 

110  CONTINUE 

C  Calculate  Seasonal  Precipitation  Patterns 

IF  (S  .eq.  1)  THEN 
CALL  Detrain 

ELSE 

CALL  Raindate 
CALL  Rainfall 

END  IF 


DO  120  Ijk  =  1,  Sumfre 

Train  =  Train  +  Rain (Ijk) 
120  CONTINUE 


C  Update  Day,  Month,  and  Equilibrium  Constants 

125  Day  =  Day  +  1 

Acumca  =  Acumca  +  Dustca 
Accums  =  Accums  +  Dusts 
IF  (Day  .gt.  Tdays)  THEN 
1  =  1  +  1 

DO  130  K  =  1,  Hori 

C02(K)  =  IC02(K)*C02Mul(I) 

130  CONTINUE 

Tdays  =  Tdays  +  Days (I) 

CALL  CONSTANTS (Ka,  Temp (Site, I) , Dhc) 

END  IF 


C  Calculate  the  Loss  of  Soil  Moisture  from  the 
C  Profile  During  the  Drying  Cycle 


140 


150 

160 

170 


Do  140  K  =  1,  Hori 

Me  =  Me  +  Moisco(K) 

CONTINUE 

IF  (Me  .gt.  Midmc)  THEN 
Devap  =  Dpotev(I) 

Me  =  0.0 

ELSE 

Soilfa  =  (Mc-Ll) / (Midmc-Ll) 

Me  =  0.0 

Devap  =  Dpotev (I) *Soilfa 

END  IF 

Tevap  =  Tevap  +  Devap 
DO  150  K  =  1,  Hori 

Availw (K) =Moisco (K) -Lowerl (K) 

IF  (Availw (K)  .gt.  Devap)  THEN 
Loss(K)  =  Devap 
Moisco(K)  =  Moisco (K) -Loss (K) 

Devap  =  0.0 
GOTO  160 

ELSE 

Loss (K)  =  Availw (K) 

Moisco (K)  =  Lowerl (K) 

Devap  =  Devap-Loss (K) 

END  IF 
CONTINUE 

IF  (Devap  .gt.  0)  Tevap  =  Tevap-Devap 

IF  ( (Rdays  .eq.  1)  .and.  (Day  .It.  365))  GOTO  125 

IF  (Day  .eq.  365)  GOTO  180 

IF  (Day  .It.  Rainda(J))  GOTO  125 

IF  (Rain(J)  .eq.  0)  THEN 
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J  =  J  +  1 
GOTO  125 

END  IF 

C  Calculate  the  Precipitation  of  CaC03  &  CaS04  During  the  Drying  Cycle 

180  DO  190  K  =  1,  Hori 

Horizo(K)  =  0.0 
Gain  (K)  =0.0 

Cas(K)  =  Ca(K)*1000.0/Moisco(K) 

S04s(K)  =  S04 (K) *1000.0/Moisco(K) 

CALL  EQUIL(K) 

Ca(K)  =  Cas(K)*Moisco{K)/1000.0 
S04(K)  =  S04s(K)*Moisco(K)/1000.0 
190  CONTINUE 

IF  ((Rainda(J)  .ne.  365)  .and.  (Day  .eg.  365))  GOTO  300 

C  Calculate  the  Gain  of  Water  During  the  Wetting  Cycle 

200  Free  =  Rain(J) 

DO  210  K  =  1,  Hori 

Availw(K)  =  Upperl (K)  -  Moisco(K) 

IF  (Availw(K)  .gt.  Rain{J))  THEN 

Moisco(K)  =  Moisco(K)  +  Rain(J) 

Gain(K)  =  Rain(J) 

Raln(J)  =  0.0 
GOTO  220 

ELSE 

Gain(K)  =  Availw(K) 

Moisco(K)  .=  Moisco(K)  +  Gain(K) 

Rain(J)  =  Rain(J)  -  Gain(K) 

END  IF 

210  CONTINUE 

220  IF  ((Rain(J)  .gt.  0)  .and.  (Cement  .eq.  0))  Leach  =  Leach+Rain ( J) 

C  Calculate  the  Rain  that  Enters  Each  Horizon 

Horizo(l)  =  Free 
DO  230  K  =  2,  Hori 

Horizo(K)  =  Horizo (K-1) -Gain (K-1) 

230  CONTINUE 

C  Calculate  the  Atmospheric  Solution  Chemistry 

Acumca  =  Acumca  +  Precca*Frec/l . 0e6 
Caa  =  Acumca/40.08 
Acumca  =  0.0 

Accums  =  Accums+Frecs*Frec/l . 0e6 
S04a  =  Accums/32 . 064 
Accums  =  0.0 

C  Mix  Horizon  Solutions  and  Reestablish  Equilibrium 
DO  240  K  =  1,  Hori 

Cas (K) = (Caa+Ca (K) ) *1000/ (Horizo (K) +Moisco (K) -Gain (K) ) 

S04s (K) = (S04a+S04 (K) ) *1000/ (Horizo (K) +Moisco (K) -Gain (K) ) 
CALL  EQUIL(K) 

Caa  =  Cas (K) * (Horizo (K) -Gain (K) ) /lOOO 
Ca(K)  =  Cas(K)*Moisco(K)/1000 
S04a  =  S04s(K)* (Horizo (K) -Gain (K))/1000 
S04(K)  =  S04s(K)*Moisco(K)/1000 
240  CONTINUE 

IF  ( (Rain (J) . gt . 0) . and. (Cement .eq. 0) )  Leacca=Leacca+Cas (K-1) * 
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X  Rain(J) *40.08/1000 

IF  ( (Rain(J) .gt.O) .and. (Cement. eq.O) )  Leachs=Leachs+S04s (K-1) * 

X  Rain(J)*32. 064/1000 

IF  (Cement . eq. 1)  CaC03 (K-1) =CaC03 (K-1) +Cas (K-1) *Rain ( J) *40 . 08/1000 
IF  (Cement. eq.l)  CaS04 (K-1) =CaS04 (K-1) +S04s (K-1) * 

X  Rain(J)*32. 064/1000 

Rain(J)  =0.0 

C  Update  Rainday 

DO  250  K  =  1,  Hori 
Gain(K)  =0.0 
Loss(K)  =0.0 
250  CONTINUE 
J  =  J  +  1 

IF  (J  .gt.  Sumfre)  THEN 
Rdays  =  1 
J  =  J  -  1 

IF  (Rainda(J)  .eq.  365)  GOTO  300 

ELSE 

IF  (Rain(J)  .eq.  0)  J=J+1 

END  IF 
GOTO  125 

C  Update  Soil  Physical  Parameters 
300  DO  310  K  =  1,  Hori 

Bd (K) =IBd (K) +CaC03 (K) *100 . 09/ (Thick (K) *40 . 08) +CaS04 (K) * 

X  172.17/ (Thick(K)*32. 064) 

IF  (Bd(K)  .gt.  2)  THEN 
Hori  =  K  -  1 
Cement  =  1 
GOTO  320 

ELSE 

Poresp(K)  =  (1.0-Bd(K)/2.65) 

IF  (Poresp(K) *Thick(K)  .gt.  Whc(K))  GOTO  305 
Upperl  (K)  =Upperl  (K)  -  (Whc  (K)  -Poresp  (K)  * 

X  Thick (K) ) 

Whc(K)  =  Upperl (K) -Lowerl (K) 

Diff  =  Whc(K)*.546 

Midmoi (K)  =  Lowerl (K)  +  Diff 

305  END  IF 

310  CONTINUE 

C  Print  Results 

320  Inte  =  Inte  +  1 

Torain  =  Torain  +  Train 
Toevap  =  Toevap  +  Tevap 
IF  (Inte  .ne.  Pint)  GOTO  350 
330  WRITE  (2,*) 

WRITE  (2,*)  ’ - 


WRITE  (2,*) 

WRITE  (2,*)  'YEAR  =',  Year 
WRITE  (2,*) 

WRITE  (2,335)  (CaC03 ( J) , J=l, Hori) 

335  FORMAT (' CAC03 (gCa/ cm2)  =',5E12.5) 
WRITE  (2,*) 

WRITE  (2,336)  (CaS04 ( J) , J=l, Hori) 

336  FORMAT (' CASO 4 (gS/ cm2)  =',5E12.5) 
WRITE  (2,*) 

WRITE  (2,337)  (Bd( J) , J=l, Hori) 
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337  FORMAT ('Bulk  Density (g/cm3)  =',5F8.4) 

WRITE  (2,*) 

WRITE  (2,345)  Leach,  Leacca,  Leachs 
345  FORMAT ('Leach (cm)  =', Ell . 4, 2X, 'Leacca (gCa/cm2)  =',E11.5,2X, 

X  'Leachs (gS/cm2)  =',E11.5) 

WRITE  (2,*) 

WRITE (2, 347)  Train, Tevap 

347  FORMAT (' Annual  Rain(cm)  =' ,F7 , 3, 5X, 'Annual  Evap(cm)  =',F7.3) 
Inte  =  0 

C  Update  Year  and  Print  Final  Profile  Concentrations 

350  Year  =  Year  +  1 

IF  (Year  .le.  Years)  GOTO  100 
WRITE  (2,*) 

WRITE  (2,*)  ' - - 


WRITE  (2,*) 

WRITE  (2,*)  'FINAL  SOIL  PROFILE  CONCENTRATIONS' 

WRITE  (2,*)  '  See  last  year  above  for  other  final  soil  properties' 

WRITE  (2,*) 

360  WRITE  (2,365) 

365  FORMAT ( ' Hori ' , 2X, 'Moisco (g/cm2) ' , 3X, ' H  activity' , 5X, 'Ca(gCa/cm2) 
x',4X, 'S04(gS/cm2) ') 

DO  370  1=1,  Hori 

Ca(I)  =  Ca(I)*40.08 
S04 (I)  =  S04 (I) *32.064 

WRITE  (2,367)  I, Moisco (I),  Hs(I),  Ca(I),  S04(I) 

367  FORMAT(lX,I2,4X,E12.5,3X,E12.5,3X,E12.5,3X,E12.5) 

370  CONTINUE 

WRITE  (2,*) 

WRITE (2, 375)  Torain,  Toevap 

375  FORMAT ('Torain (cm)  =' , 1X,E15. 6, 5X, ' Toevap (cm)  =',1X,E15.6) 

WRITE  (2,*) 

WRITE (2, 376)  Acumca, Accums 

376  FORMAT ( 'Residual  Dust  Ca{gCa/cm2)  =', Ell. 4, 5X, 'Residual  Dust  S04 
x(gS/cm2)  =',E11.4) 

END 

C - 

C - 

SUBROUTINE  Detrain 

C  Subroutine  for  Deterministic  Rain  Model.  The  Dimension  statement, 

C  DATA,  and  Sumfre  must  be  Changed  for  each  Specific  Site. 

DIMENSION  Raindy (45) ,  Raina(45) 

INTEGER  Sumfre,  Site 

COMMON  /  Precipitation  /  Rain (200),  Rainda(200),  Sumfre,  Site, 

X  Delrai,  Numb,  Rnum,  Inter (7, 3, 11) ,  Freq(7, 3, 11) ,  R(7,3,ll), 

X  D(7,3,ll) 

DATA  Raindy  /  9,10,11,19,21,29,31,39,44,45, 

X  46,47,49,52,66,69,70,78,85,86, 

X  120,172,180,181,187,193,200,204,206,208, 

X  210,214,225,226,235,236,248,250,267, 

X  268,269,288,289,340,346  / 

DATA  Raina  /  . 03, . 05, . 10, . 61, . 61, . 15, . 05, 2 . 18, 2 . 11, . 99, 

X  1.14, .05, .05, .05, .10, .48, .36, .38,1.47, .30, 

X  .20, .20, .03, .36, .25,1.60, .84, .03, .38, .10, 

X  .18, .08,1.75,1.09, .18,1.47, .74,2.74, .99, 

X  1.98, .91, .05, .51, .10, .38  / 

Sumfre  =45 
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DO  440  Ijk  =  1,  Sumfre 

Rainda(Ijk)  =  Raindy(Ijk) 

Rain(Ijk)  =  Raina(Ijk) 

440  CONTINUE 

RNumb  =  Delrai* (Sumfre) 

Numb  =  NINT (RNumb) 

DO  450  I  =  1,  Numb 
RI  =  I 

RNumb  =  Numb 
Sum  =  Sumfre 
Elim  =  RI*Sum/RNumb 
Elim  =  INT(Elim) 

Rain (Elim)  =  0 
450  CONTINUE 
RETURN 
END 

C - 

C - 

SUBROUTINE  Raindate 

C  Calculate  rain  dates  using  a  stochastic  model 
C  Present  model  developed  for  seven  southwestern  sites 
C  Frequency  distributions  need  to  be  changed  for  other  sites 

REAL  Newday (100) ,  S(10),  Intday 
REAL* 4  Rnum,  rand 

INTEGER  Site,  Season,  X,  Y,  Sumfre 

COMMON  /  Precipitation  /  Rain (200) , Rainda (200) , Sumfre, Site, 

X  Delrai,  Numb,  Rnum,  Inter  (7,  3, 11) ,  Freq  (7,  3, 11) ,  R  (7, 3, 11) , 

X  D(7,3,ll) 

C  Calculate  Raindates 

J  =  1 
Day  =  0 

485  CALL  SEASONS (Site, Season, Day, Numb, New j) 

DO  490  1=1,  10 

S (I) = (Inter (Site, Season, I+l) -Inter (Site, Season, I) ) / 

X  (Freq (Site, Season, I+l) -Freq (Site, Season, I)  ) 

490  CONTINUE 

Rnum  =  rand(O) 

DO  500  1=2,  11 

IF  (Rnum  .It.  Freq(Site, Season, I) )  GOTO  510 
500  CONTINUE 

510  lntday=lnter (Site, Season, I-l) +S (I-l) * (Rnum-Freq (Site, Season, I-l) ) 
x+1.000 

Day  =  Day+Intday 

IF  (Day  .gt.  365)  GOTO  515 

Rainda (J)  =  Day 

J  =  J+1 

GOTO  485 

515  J  =  J-1 

IF  (Newj  .eq.  0)  GOTO  555 

C  Add  new  raindates  to  reproduce  long-term  mean  (if  necessary) 

C  and  sort  dates  numerically 

DO  520  1=1,  Newj 

Rnum  =  rand(O) 

Newday (I)  =  Rnum*364+1 . 000 


26 


4/18/94  10:09 


Macintosh  HD : MPW ; CALGYP . f 


520  CONTINUE 

J  =  J+Newj 
Ja  =  J+l-Newj 
1  =  1 

DO  530  K  =  Ja,  J 

Raincia(K)  =  Newday(I) 

I  =  I+l 

530  CONTINUE 
Jb  =  J-1 
DO  550  X  =  1,  J 

DO  540  Y  =  1,  Jb 

IF  (Rainda(Y)  .le.  Rainda(Y+l))  GOTO  540 
Temper  =  Rainda(Y) 

Rainda(Y)  =  Rainda(Y+l) 

Rainda(y+1)  =  Temper 
540  CONTINUE 

550  CONTINUE 

C  Eliminate  Duplicate  Raindates 

555  Dup  =  0 
1  =  2 

560  N  =  I+l 

IF  (INT(Rainda(I) )  .eq.  INT (Rainda (I-l) ) )  THEN 
Dup  =  Dup+1 

IF  (N  .gt.  J)  GOTO  575 
DO  570  K  =  N,  J 

Rainda (K-1)  =  Rainda (K) 


570 

CONTINUE 

575 

J  =  J-1 

END  IF 
I  =  I  +  l 

IF  (I  .le.  J)  GOTO  560 
IF  (Dup  .gt.  0)  GOTO  555 

C  Return  to  CALGYP  with  New  Raindates 

Sumfre  =  J 

DO  580  1=1,  Sumfre 

Rainda (I)  =  INT (Rainda (I) ) 

580  CONTINUE 
RETURN 
END 

C - 

C - 

SUBROUTINE  Rainfall 

C  Calculate  the  storm  rainfall  amounts  using  a  stochastic  model 
C  Present  model  developed  for  seven  southwestern  sites 
C  Frequency  distributions  need  to  be  changed  for  other  sites. 

REAL  S(10) 

REAL* 4  Rnum,  rand 

INTEGER  Season,  Site,  Sumfre 

COMMON  /  Precipitation  /  Rain (200) , Rainda (200) , Sumfre, Site, 

X  Delrai,Numb, Rnum, Inter (7, 3, 11) , Freq(7, 3,11),R(7,3,11), 

X  D(7,3,ll) 

C  Calculate  Rainfall  Amounts 

J  =  Sumfre 


27 


4/18/94  10:09 


Macintosh  HD : MPW : CALGYP . f 


K  =  1 

600  CALL  SEASONS (Site,  Season,  Rainda(K),  Numb,  Newj) 

DO  610  1=1,  10 

S (I ) = (R (Site, Season, I+l ) -R (Site, Season, I) ) / 

X  (D  (Site,  Season,  I+l)  -D  (Site,  Season,  I) ) 

610  CONTINUE 

Rnum  =  rand(O) 

DO  620  1=2,  11 

IF  (Rnum  .It.  D (Site, Season, I) )  GOTO  630 
620  CONTINUE 

630  Rain(K)  =  R(Site, Season, I-l) +S (I-l) * (Rnum-D (Site, Season, I-l) ) 
K=K+1 

IF  (K  .le.  J)  GOTO  600 

C  Set  subset  of  rainfall  amounts  to  zero  to  reproduce  long-term 
C  mean  (if  necessary)  and  to  produce  drier  climates. 

RNumb  =  Numb  +  Delrai* (J-Numb) 

Numb  =  NINT (RNumb) 

Do  640  1=1,  Numb 
RI  =  I 
RJ  =  J 
RNumb  =  Numb 
Elim  =  RI*RJ/RNumb 
Elim  =  INT(Elim) 

Rain (Elim) =0.0 
640  CONTINUE 

C  Return  to  CALGYP  with  New  Rainfall  Amounts 

RETURN 

END 

C - 

C - 

SUBROUTINE  SEASONS (Site, Season, Day, Numb, New j ) 

C  This  subroutine  calculates  the  proper  season  for  a  given  day  and  site 

INTEGER  Site,  Season 

IF  (Site  .eq.  1)  THEN 
Numb=7 
New j=0 

IF  ((Day  .It.  152)  .OR.  (Day  .gt.  304))  THEN 
Season  =  1 

ELSE 

Season  =  2 

END  IF 

END  IF 

IF  (Site  .eq.  2)  THEN 
Newj  =  7 
Numb=0 

IF  ((Day  .It.  182)  .OR.  (Day  .gt.  304))  THEN 
Season  =  1 

ELSE 

Season  =  2 

END  IF 

END  IF 

IF  (Site  .eq.  3)  THEN 
Newj  =  6 
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Numb=0 

IF  {(Day  .It.  121)  .OR.  (Day  .gt.  273))  THEN 
Season  =  1 

ELSE 

Season  =  2 

END  IF 

END  IF 

IF  (Site  .eq.  4)  THEN 
Nurtib=3 
Newj=0 

IF  ((Day  .It.  121)  .OR.  (Day  .gt.  304))  THEN 
Season  =  1 

ELSE 

Season  =  2 

END  IF 

END  IF 

IF  (Site  .eq.  5)  THEN 
Numb=l 
New j=0 

IF  ((Day  .It.  91)  .OR.  (Day  .gt.  212))  THEN 
Season  =  1 

ELSE 

Season  =  2 

END  IF 

END  IF 

IF  (Site  .eq.  6)  THEN 
Numb=4 
New j=0 

IF  ((Day  .It.  91)  .OR.  (Day  .gt.  273))  THEN 
Season  =  1 

ELSE 

IF  (Day  .It.  182)  THEN 
Season  =  2 

ELSE 

Season  =  3 

END  IF 

END  IF 

END  IF 


IF  (Site  .eq.  7)  THEN 
Numb=9 
New j=0 

IF  {(Day  .It.  91)  .OR.  (Day  .gt.  273))  THEN 
Season  =  1 

ELSE 

IF  (Day  .It.  182)  THEN 
Season  =  2 

ELSE 

Season  =  3 

END  IF 

END  IF 

END  IF 


RETURN 

END 


C 


SUBROUTINE  CONSTANTS (Ka, Temp, Dhc) 
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C  Calculate  the  Equilibrium  Constants  as  Functions  of  Temperature 

DOUBLE  PRECISION  Ka{6) 

DIMENSION  Pk(6) 

Pk(l)  =  l,14+.0131*Temp 
Pk(2)  =  6. 54-. 0071*Temp 
Pk(3)  =  10.59-.0102*Temp 
PkU)  =  7.95+.0125*Temp 
Pk(5)  =  4.62+.0006*Temp 
Pk(6)  =  2.23+. 0019*Temp 
DO  480  K  =  1,  6 

Ka(K)  =  EXP(-2.3026*Pk(K)) 

480  CONTINUE 

Dhc=0. 4918+6. 60 98e-4*Temp+5.0231e-6*Temp**2 

RETURN 

END 

C - 

C - 

SUBROUTINE  ION (lonstr , Acl , Ac2, Dhc) 

C  Calculate  Uni-  and  Di-valent  Activity  Coefficients  Using 
C  the  Davies  Equation 

REAL  lonstr 

Factor  =  SORT (lonstr )/(1.0+SQRT (lonstr ))-.3*Ionstr 
Acl  =  EXP(-2.3026*Dhc*Factor) 

Ac2  =  EXP(-4*2.3026*Dhc*Factor) 

RETURN 

END 

C - 

C - 

SUBROUTINE  Hact  (C02,  Hs,  S04,  Acl,  Ac2,  Ka) 

C  This  is  a  subroutine  to  calculate  the  H  activity  for  a 
C  system  in  equilibrium  with  CaC03  and  CaS04. 

DOUBLE  PRECISION  Ka (6) , Kl, K2, K3, F, Df , H, Dx, Hs, S04 
DIMENSION  Frac(8) 

H=Hs 

DATA  Frac/.8,1.2, . 5, 1 . 5, . 2, 5 . 0, . 1, 10 . 0/ 

1=1 

Kl=2 . 0*Ka (4) / (Ka (3) *Ka (2) *Ka (1) ) 

K2=2.0*Ka(l)*Ka(2)*Ka(3) 

K3=Ka(l) *Ka(2) 

375  F=K1*H**2/ (C02*Ac2) -K2*C02/ (H**2*Ac2) -K3*C02/ (H*Acl) -2 . 0*SO4 

Df=2 . 0*K1*H/ (C02*Ac2) +K3*C02/ (H**2*Acl) +2*K2*C02/ (H**3*Ac2) 
F=-1.0*F 
Dx=F/Df 
H=H+Dx 

Pcen=ABS (Dx/H) *100.0 
IF  (Pcen  .It.  1)  GOTO  380 

IF  ((H  .gt.  l.Oe-5)  .or.  (H  .It.  l.Oe-10))  THEN 
H=Hs*Frac(I) 

1=1+1 

IF  (I  .gt.  8)  THEN 

PRINT*, 'The  hydrogen  subroutine  will  not  converge' 
GOTO  390 

END  IF 

END  IF 


30 


4/18/94  10:09 


Macintosh  HD : MPW : CALGYP . f 


GOTO  375 
380  Hs=H 

RETURN 
390  END 

C - 

C - 

SUBROUTINE  EQUIL (K) 

C  This  is  a  subroutine  to  determine  equilibrium  concentrations  of 
C  solutes  in  saturated  CaC03  and  CaS04  solutions 

REAL  I 

DOUBLE  PRECISION  A, B, C, X, Caion, S04ion, Cas, S04s, CaSO40, SysCa, 
xCaC03,Pk4,Pk5,SysS04,CaS04,Pk6,Test,HC03,C03,Ka,Hs,Moisco, 
xHorizo, Gain 

COMMON  /  Chemistry  /  SO4s(10),  CaCO3(10),  CaSO4(10),  Cas(10),Dhc, 
X  CO2(10),  Hs(lO),  Ka(6),  Moisco(IO),  Horizo(lO),  Gain(lO) 

C  Initiate  Variables 


I  =  3.0*Cas(K) 
Caion  =  Cas(K) 
Test  =  Caion 


IF  ((S04s(K)  .EQ.  0.0)  .AND. 

(CaS04 (K) 

.EQ. 

0.0))  THEN 

SO41on=0.0 

CaS040=0.0 

END  IF 

460 

CALL  ION (I,  Acl,  Ac2,Dhc) 

IF  ((S04s(K)  .EQ.  0.0)  .AND. 

(CaS04 (K) 

.EQ. 

0.0) )  GO  TO  465 

C  Calculate  CaS04  lonpair 

Pk6  =  Ka(6) /Ac2**2 
A  =  1.0 

B  =  -l,0*(Cas(K)+SO4s(K)+Pk6) 

C  =  Cas(K)*S04s(K) 

X=  (-B-SQRT(B**2-4.0*A*C) )/2.0*A 
Caion  =  Cas (K)  -  X 
S04ion  =  S04s(K)  -  X 
CaSO40  =  X 


C  Calculate  CaC03  Precipitation-Dissolution 

465  IF  (CaC03(K)  .gt.  0.0)  CALL  Hact {C02 (K) , Hs (K) , S04ion, Acl, Ac2, Ka) 
SysCa=Caion+CaC03 (K) *1000.0/(40. 08* (Moisco (K) +Horizo (K) - 
X  Gain (K) ) ) 

Pk4= (Hs (K) **2*Ka (4) ) / (Ac2*Ka (1) *Ka (2) *Ka (3) *C02 (K) ) 

X  =  SysCa-Pk4 

IF  (X  .It.  0.0)  X=0.0 

Caion  =  SysCa-X 

CaC03 (K) = (X) *40 . 08* (Moisco (K) +Horizo (K) -Gain (K) ) /lOOO . 0 
IF  ((S04s(K)  .EQ.  0.0)  .AND.  (CaS04 (K)  .EQ.  0.0))  GO  TO  467 

C  Calculate  CaS04  Precipitation-Dissolution 

SysCa=Caion+CaS04 (K) *1000 . 0/ (32 . 064* (Moisco (K) +Horizo (K) - 
X  Gain (K) ) ) 

SysS04=S04ion+CaS04 (K) *1000 . 0/ (32 . 064* (Moisco (K) + 

X  Horizo (K) -Gain (K) ) ) 

Pk5  =  Ka(5)/Ac2**2 
A  =  1.0 
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B  =  -1.0* (SysCa+SysS04) 

C  =  SysCa*SysS04-Pk5 
X=  (-B-SQRT (B**2-4*A*C) ) / (2.0*A) 

IF  (X  .It.  0.0)  X  =  0.0 
Caion  =  SysCa-X 
S04ion  =  SysS04-X 

CaS04 (K) = (X) *32 . 064* (Moisco (K) +Horizo (K) -Gain (K) ) /lOOO . 0 

C  Update  Concentrations  and  Test  for  Convergence 

467  Cas(K)  =  Caion+CaSO40 

S04s (K)  =  SO4ion+CaSO40 

HC03  =  Ka(l)*Ka(2)*C02(K)/(Hs(K)*Acl) 

C03=Ka  (1)  *Ka  (2)  *Ka  (3)  *C02  (K)  /  (Hs  (X)  **2*Ac2) 

I  =  .5* (4.0* (Caion+S04ion+C03)+HC03) 

Diff  =  ABS((Test-Caion) /Caion) *100.0 
IK  (Diff  .gt.  1.0)  THEN 
Test  =  Caion 
GOTO  460 

END  IF 

470  RETURN 
END 

C - 
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PROGRAM  RAINMODULE 

This  program  calculates  raindates  and  rainfall  amounts  using  a 
stochastic  model  for  seven  southwestern  sites. 

INTEGER  Site,  Year,  Sum,  Sumfre,  Totsum,  Seed 
REAL  *4  Rnum,  rand 
CHARACTER*50  Title 

COMMON  Rain (200) ,  Rainda(200),  Suipfre,  Site,  Delrai,  Numb, Rnum, 
X  Inter  (7,3,11) ,  Freq (7, 3, 11) ,  R(7,3,ll),  D(7,3,ll) 

C  Read  in  the  program  data 

PRINT*,  'Select  Site:' 

PRINT*,  '  1  =  El  Paso,  TX  (21.6cm)' 

PRINT*,  '  2  =  Albuq.,  NM  (21.1cm)' 

PRINT*,  '  3  =  Clayton,  NM  (37.8cm)' 

PRINT*,  '  4  =  Roswell,  NM  {31.6cm)' 

PRINT*,  '  5  =  Yuma,  AZ  (8.5cm)' 

PRINT*,  '  6  =  Phoenix,  AZ  (18.9cm)' 

PRINT*,  '  7  =  Tucson,  AZ  (28.4cm)' 

READ*,  Site 

C  Read  in  basic  model  run  parameters 

PRINT*,  'Title?' 

READ*,  Title 

PRINT*,  'Fractional  Change  in  Rain  for  Drier  Climate.' 

PRINT*,  '  Default  =  0.00,  no  change.' 

READ*,  Delrai 

PRINT*,  'Seed  for  Random  Number  Generator?' 

READ*,  Seed 
Rnum  =  rand (Seed) 

PRINT*,'  Rain(cm)  Rain  Events' 

C  Read  in  Stochastic  Rain  Model  Parameters 


DATA 

(((D(I 

,  J,K) 

,K=1, 

11), 

J=1,2),I=1,7)  / 

X 

.000, 

.571, 

.745, 

.827 

,  .837, .888, .918 

,  .939 

,  .959, .990, 

1. 

000, 

EP, 

WINT 

X 

.000, 

.526, 

.691, 

.773 

,  .835, .866,  .897 

,  .938 

,  .948, .979, 

1. 

000, 

EP, 

SUM 

X 

.000, 

.  635, 

.832, 

.920 

,  .964, .971, .986 

,  .993 

,1.000,1.05 

,1 

-.10, 

AL, 

WINT 

X 

.000, 

.577, 

.740, 

.827 

,  .904,  .952, .962 

,  .981 

,  .991,1.00, 

1. 

05, 

AL, 

SUM 

X 

.000, 

.586, 

.811, 

.928 

,  .955, .982, .991 

,1.00 

,1.05,1.10, 

1. 

15, 

CL, 

WINT 

X 

.000, 

.371, 

.596, 

.728 

,  .794, .860, .880 

,  .933 

,  .966, .986, 

1. 

000, 

CL, 

SUM 

X 

.000, 

.570, 

.756, 

.884 

,  .907, .942, .965 

,1.00 

,1.05,1.10, 

1. 

15, 

RS, 

WINT 

X 

.000, 

.405, 

.628, 

.719 

,  .769, .827, .868 

,  .901 

,  .975, .983, 

1. 

000, 

RS, 

SUM 

X 

.000, 

.512, 

.780, 

.865 

,  .889, .938, .962 

,  .986 

,1.00,1.05, 

1. 

10, 

YU, 

WINT 

X 

.000, 

.625, 

.875, 

1.000,1.05,1.10,1.15,1.20,1.25,1.30 

-.35, 

YU, 

SUM 

X 

.000, 

.394, 

.567, 

.711 

, .817, .875, .923 

,  .952 

, .971, .990, 

1. 

,000, 

PH, 

WINT 

X 

.000, 

.600, 

.800, 

.933 

,1.000,1.1,1.2, 

1.3,1 

.4, 1.5, 1.6, 

PH, 

SPR 

X 

.000, 

.473, 

.661, 

.732 

,  .803, .830, .875 

,  .902 

, .947, .983, 

1. 

000, 

TU, 

WINT 

X 

.000, 

.550, 

.850, 

.950 

,1.000,2.0,3.0, 

4.0,5 

.0,6. 0,7.0 

/ 

TU, 

SPR 

DATA 

(  (D(I, 

3,K)  , 

K=l,ll)  ,I 

=6,7)  / 

X 

.000, 

.642, 

.755, 

.793 

, .868, .925, .963 

,  .982 

,1.000,1.5, 

2. 

.0, 

PH, 

SUM 

X 

.000, 

.426, 

.539, 

.617 

, .730, .765, .843 

,  .878 

, .913, .991, 

1. 

.000  / 

TU, 

SUM 

DATA 

( ( (R(I 

,  J,K) 

,K=1, 

11), 

J=1,2),I=1,7)  / 

X 

0.00, 

.25,  . 

50, .75,1. 

00,1.25,1.50,1. 

75,2. 

00, 3.00, 3.07, 

EP, 

WINT 

X 

0.00, 

.25,  . 

50, .75,1. 

00,1.25,1.50,1. 

75,2. 

00,4.00,5.59, 

EP, 

SUM 

X 

0.00, 

.25,  . 

50, .75,1. 

00,1.25,1.50,2. 

00,2. 

64,3.00,5.00, 

AL, 

WINT 

X 

0.00, 

.25,  . 

50, .75,1. 

00,1.25,1.50,2. 

00,3. 

00,4.45,10. 

0, 

AL, 

SUM 

X 

0.00, 

.25,  . 

50,  .7 

'5,1. 

00,1.25,2.00,4. 

83,5. 

00,10.0,15. 

0, 

CL, 

WINT 
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X  0.00, .25, .50, .75,1.00,1,50,2.00,3.00,4.00,5.00,5.36,  CL, SUM 

X  0.00, .25, .50, .75,1.00,1,25,1.50,1,75,5.00,10.0,15.0,  RS,WINT 

X  0.00, .25, .50, .75,1.00,1.50,2.00,3.00,5.00,7.00,10.72,  RS,SUM 

X  0.00, .25, .50, .75,1.00,1.25,1,50,3.00,4,45,10.0,16.0,  YU,WINT 

X  0.00, .25, .50, .58,1.00,1.25,1.50,1.75,2.00,3.00,3.07,  YU, SUM 

X  0.00, .25, .50, .75,1.00,1.25,1.50,1,75,2.00,3.00,5.03,  PH,WINT 

X  0.00, .25, .50, .75, .97,1.25,1.50,1.75,2.00,3.00,3.07,  PH,SPR 

X  0.00, .25, .50, .75,1.0,1.25,1.5,1.75,2.0,3.0,7.52,  TU,WINT 

X  0.00, .25, .50, .75,1.02,2.0,3.0,4.0,5.0,6.0,7.0  /  TU,SPR 

DATA  ((R{I,3,K),K=1,11),I=6,7)  / 

X  0.00, .25, .50, .75,1.00,1.50,2.00,3.50,4.22,5.00,10.0,  PH, SUM 

X  0.00, .25, .50, .75,1.00,1.25,1.50,1.75,2.00,3.00,3.23  /  TU,SUM 


DATA  (  (  (Freqd,  J,K)  ,K=1,11),  J=1,2),I=1,7)  / 

X  0.00, .388, .449, .510, .531, .643, .765, .796, .888, .969,1.000, 

X  0.00, .402, .454, .567, .680, .750, .853, .902, .923, .974,1.000, 

X  0.00, .328, .447, .522, . 589, . 664, . 768, .843, .910, .977,1.000, 

X  0.00, .396, .481, .594, .679, .783, .840, .963, .982, .991,1.000, 

X  0.00, .291, .409, .427, .472, .617, .762, .844, .944, .989,1.000, 

X  0.00, .404, .530, .623, .716, .848, .901, .961, .987, .994,1.000, 

X  0.00, .430, .488, .569, .616, .744, .825, .918, .988,1.000,1.05, 

X  0.00, .380, .430, .521, .554, .703, .802, .918, .984, .992,1,000, 

X  0.00, .316, .417, .430, ,493, .582, .671, .785, .950, .988,1.000, 

X  0.00, .111, .333, .777,1.000,1.05,1.1,1.2,1.3,1.4,1,5, 

X  0.00, .476, .515, ,583, .641, .719, .787, .865, ,923, .972,1.000, 

X  0.00, .250, .313, ,438, .501, .564, .689, .814,  .939,1.00,1.05, 

X  0.00, .447, .552, .596, .631, .745, .815, .885, .955, .973,1.000, 

X  0.00, .222, .278, .334, .390, .501, .612, .779, .890,  .946,1.000  / 

DATA  ((Freq(l,3,K),K=l,ll),I=6,7)  / 

X  0.00, .212, .347, .482, .559, .713, .790, .886, .963,1.00,1.05,  PH, SUM 

X  0.00, .491, .640, .728, .781, .886, .965, .991,1.00,1.1,1.2  /  TU, SUM 

DATA  (((Interd,  J,K),K=1,11),J=1,2),I=1,7)  / 

X  0,1,2,3,4,7,10,15,20,30,57, 

X  0,  1,  2,  3,  4,7,10,15,20,30,  66, 

X  0,1,2,3,4,7,10,15,20,30,40, 

X  0,1,2,3,4,7,10,15,20,30,32, 

X  0,1,2,3,4,7,10,15,25,40,50, 

X  0,1,2,3,4,7,10,13,15,20,21, 

X  0,1,2,3,7,10,15,25,40,64,100, 

X  0,1,2,3,4,7,10,15,25,40,48, 

X  0,1,2,3,4,7,10,15,40,70,111, 

X  0,20,40, 60,121,150,160,170,180,190,200, 

X  0,1,2,3,4,7,10,15,25,40,52, 

X  0,1,5,10,15,25,40,50,60,77,100, 

X  0,1,2,3,4,7,10,15,25,40,57, 

X  0,1,3,5,7,10,25,30,40,60,78  / 

DATA  ( (Inter (1,3, K) ,K=1, 11) ,1=6,7)  / 

X  0,1,2,3,4,7,10,15,20,24,50,  PH, SUM 

X  0,1,2,3,4,7,10,15,17,40,50  /  TU, SUM 

Inte  =  1 
Torain  =  0.0 
Totsum  =  0.0 

10  Year  =  1 

Cenrai  =  0.0 
Sum  =  0 

20  Train  =  0.0 

C  Calculate  Annual  Rainfall  Pattern 


EP,WINT 

EP,SUM 

AL,WINT 

AL, SUM 

CL,WINT 

CL, SUM 

RS,WINT 

RS,SUM 

YU,WINT 

YU, SUM 

PH,WINT 

PH, SPR 

TU,WINT 

TU, SPR 


EP,WINT 
EP,SUM 
AL,WINT 
AL, SUM 
CL,WINT 
CL, SUM 
RS,WINT 
RS,SUM 
YU, HINT 
YU, SUM 
PH,WINT 
PH, SPR 
TU,WINT 
TU, SPR 
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CALL  Raindate 
CALL  Rainfall 

DO  30  Ijk  =  1,  Sumfre 

Train  =  Train  +  Rain (Ijk) 

30  CONTINUE 

Sumfre  =  Sumfre  -  Numb 
Cenrai  =  Cenrai  +  Train 
Torain  =  Torain  +  Train 
Sum  =  Sum  +  Sumfre 
Totsum  =  Totsum  +  Sumfre 

C  Update  Annual  Parameters 

Year  =  Year  +  1 

IF  (Year  .LE.  100)  GOTO  20 

50  PRINT  60,  Cenrai,  Sum 

60  FORMAT  (IX, '100  yrs  =  ',  F8.1,  3X,  16) 

Inte  =  Inte  +  1 

IF  (Inte  .LE.  10)  GOTO  10 

70  PRINT  80,  Torain,  Totsum 

80  FORMAT  ('1000  yrs  =  ',  F8.1,  3X,  16) 

PAUSE 

END 

C - 

C - 

SUBROUTINE  Raindate 

C  Calculate  rain  dates  using  a  stochastic  model 
C  Present  model  developed  for  seven  southwestern  sites 
C  Frequency  distributions  need  to  be  changed  for  other  sites 

REAL  Newday (100) ,  S(10),  Intday 
REAL* 4  Rnum,  rand 

INTEGER  Site,  Season,  X,  Y,  Sumfre 
COMMON  Rain (200) ,Rainda(200) , Sumfre, Site, 

X  Delrai,  Numb, Rnum, Inter (7, 3, 11) , Freq (7, 3, 11) , R(7, 3, 11) , 

X  D(7,3,ll) 

C  Calculate  Raindates 

J  =  1 
Day  =  0 

485  CALL  SEASONS (Site, Season, Day, Numb, Newj) 

DO  490  1=1,  10 

S (I) = (Inter (Site, Season, I+l) -Inter (Site, Season, I) ) / 

X  (Freq(Site, Season, I+l)-Freq(Site, Season, I) ) 

490  CONTINUE 

Rnum  =  rand(O) 

DO  500  1=2,  11 

IF  (Rnum  .It.  Freq(Site, Season, I) )  GOTO  510 
500  CONTINUE 

510  Intday=Inter (Site, Season, I-l) +S (I-l) * (Rnum-Freq (Site, Season, I-l) ) 
x+1.000 

Day  =  Day+Intday 
IF  (Day  .gt.  365)  GOTO  515 
Rainda(J)  =  Day 
J  =  J+1 
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GOTO  485 
515  J  =  J-1 


IF  (Newj  .eq.  0)  GOTO  555 

C  Add  new  raindates  to  reproduce  long-term  mean  (if  necessary) 
C  and  sort  dates  numerically 


DO  520  1=1,  Newj 

Rnum  =  rand(O) 

Newday(I)  =  Rnum*364+1 . 000 
520  CONTINUE 

J  =  J+Newj 
Ja  =  J+l-Newj 
1  =  1 

DO  530  K  =  Ja,  J 

Rainda(K)  =  Newday(I) 

I  =  I+l 

530  CONTINUE 
Jb  =  J-1 
DO  550  X  =  1,  J 

DO  540  Y  =  1,  Jb 

IF  (Rainda(Y)  .le.  Rainda(Y+l))  GOTO  540 
Temper  =  Rainda(Y) 

Rainda(Y)  =  Rainda(Y+l) 

Rainda(Y+l)  =  Temper 
540  CONTINUE 

550  CONTINUE 


C  Eliminate  Duplicate  Raindates 


555  Dup  =  0 
1  =  2 

560  N  =  I+l 

IF  (INT(Rainda(I) )  .eq.  INT (Rainda (I-l) ) )  THEN 
Dup  =  Dup+1 

IF  (N  .gt.  J)  GOTO  575 
DO  570  K  =  N,  J 

Rainda (K-1)  =  Rainda (K) 

570  CONTINUE 


575  J  =  J-1 

END  IF 
I  =  I+l 

IF  (I  .le.  J)  GOTO  560 
IF  (Dup  .gt.  0)  GOTO  555 


C  Return  to  Rainmodule  with  New  Raindates 


Sumfre  =  J 

DO  580  1=1,  Sumfre 

Rainda (I)  =  INT (Rainda (I) ) 
580  CONTINUE 
RETURN 
END 


SUBROUTINE  Rainfall 

C  Calculate  the  storm  rainfall  amounts  using  a  stochastic  model 
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C  Present  model  developed  for  seven  southwestern  sites 
C  Frequency  distributions  need  to  be  changed  for  other  sites. 

REAL  S(10) 

REAL*4  Rnum,  rand 

INTEGER  Season,  Site,  Sumfre 

COMMON  Rain (200) ,Rainda (200) , Sumfre, Site, 

X  Delrai,Numb, Rnum, Inter (7, 3, 11) , Freq(7, 3,11),R(7,3,11), 

X  D(7,3,ll) 

C  Calculate  Rainfall  Amounts 

J  =  Sumfre 
K  =  1 

600  CALL  SEASONS (Site,  Season,  Rainda(K),  Numb,  Newj) 

DO  610  1=1,  10 

S (I) = (R (Site, Season, I+l) -R(Site, Season, I) ) / 

X  (D (Site, Season, I+l) -D (Site, Season, I) ) 

610  CONTINUE 

Rnum  =  rand(O) 

DO  620  1=2,  11 

IF  (Rnum  .It.  D (Site, Season, I) )  GOTO  630 
620  CONTINUE 

630  Rain(K)  =  R (Site, Season, I-l) +S (I-l) * (Rnum-D (Site, Season, I-l) ) 
K=K+1 

IF  (K  .le.  J)  GOTO  600 

C  Set  subset  of  rainfall  amounts  to  zero  to  reproduce  long-term 
C  mean  (if  necessary)  and  to  produce  drier  climates. 

RNumb  =  Numb  +  Delrai* ( J-Numb) 

Numb  =  NINT (RNumb) 

Do  640  1=1,  Numb 
RI  =  I 
RJ  =  J 
RNumb  =  Numb 
Elim  =  RI*RJ/RNumb 
Elim  =  INT(Elim) 

Rain (Elim) =0 . 0 
640  CONTINUE 

C  Return  to  Rainmodule  with  New  Rainfall  Amounts 

RETURN 

END 


C - 

C - 

SUBROUTINE  SEASONS (Site, Season, Day, Numb, New j ) 

C  This  subroutine  calculates  the  proper  season  for  a  given  day  and  site 

INTEGER  site, Season 

IF  (Site  .eq.  1)  THEN 
Numb=7 
New j=0 

IF  ((Day  .It.  152)  .OR.  (Day  .gt.  304))  THEN 
Season  =  1 

ELSE 

Season  =  2 

END  IF 
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END  IF 

IF  (Site  .eq.  2)  THEN 
Newj  =  7 
Numb=0 

IF  ((Day  .It.  182)  .OR.  (Day  .gt.  304)) 
Season  =  1 

ELSE 

Season  =  2 

END  IF 

END  IF 

IF  (Site  .eq.  3)  THEN 
Newj  =  6 
Numb=0 

IF  ((Day  .It.  121)  .OR.  (Day  .gt.  273)) 
Season  =  1 

ELSE 

Season  =  2 

END  IF 

END  IF 

IF  (Site  .eq.  4)  THEN 
Niimb=3 
New j=0 

IF  ((Day  .It.  121)  .OR.  (Day  .gt.  304)) 
Season  =  1 

ELSE 

Season  =  2 

END  IF 

END  IF 

IF  (Site  .eq.  5)  THEN 
Nuinb=l 
New j=0 

IF  ((Day  .It.  91)  .OR.  (Day  .gt.  212)) 
Season  =  1 

ELSE 

Season  =  2 

END  IF 

END  IF 

IF  (Site  .eq.  6)  THEN 
Numb=4 
New j=0 

IF  ((Day  .It.  91)  .OR.  (Day  .gt.  273)) 
Season  =  1 

ELSE 

IF  (Day  .It.  182)  THEN 
Season  =  2 

ELSE 

Season  =  3 

END  IF 

END  IF 

END  IF 

IF  (Site  .eq.  7)  THEN 
Nuinb=9 
New j=0 

IF  ((Day  .It.  91)  .OR.  (Day  .gt.  273)) 
Season  =  1 

ELSE 

IF  (Day  .It.  182)  THEN 


THEN 


THEN 


THEN 


THEN 


THEN 


THEN 
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Macintosh  HD:MPW:RainMd. f 


Season  =  2 

ELSE 

Season  =  3 

END  IF 

END  IF 

END  IF 
RETURN 
END 
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